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IMPLEMENTATION  OF  WETTING  AND  DRYING  IN  NCOM: 
DESCRIPTION  AND  VALIDATION  TEST  REPORT 


1.  INTRODUCTION 

Wetting  and  Drying  (WAD)  refers  to  the  wetting  or  flooding  of  grid  cells  that  were  previously 
dry  and  the  drying  of  grid  cells  that  were  previously  wet  (i.e. ,  covered  by  water)  during  the  temporal 
integration  of  an  ocean  or  fluid-flow  model.  The  Navy  Coastal  Ocean  Model  (NCOM)  (Martin  2000, 
Morey  et  al.  2003,  Barron  et  al.  2004)  does  not  currently  have  a  WAD  capability;  hence,  grid  cells 
cannot  switch  between  being  wet  and  dry,  i.e.,  grid  cells  that  are  initially  dry  must  stay  dry  and 
those  that  are  initially  wet  must  stay  wet.  There  is  no  mechanism  by  which  the  WAD  “front” 
between  dry  grid  cells  and  wet  grid  cells  can  advance  onto  the  dry  cells  when  the  sea-surface  height 
(SSH)  at  the  coast  rises.  Similarly,  there  is  no  mechanism  by  which  the  WAD  front  between  the 
dry  grid  cells  and  the  wet  grid  cells  can  retreat  from  the  previously  wet  grid  cells  when  the  SSH  at 
the  coast  falls. 

The  importance  of  WAD  in  a  simulation  depends  both  on  the  extent  of  the  WAD  areas  in  the 
domain  covered  by  the  ocean  model  and  the  grid  resolution  being  used. 

The  rise  of  sea  level  at  wet  grid  cells  near  the  coast  in  NCOM  does  not  cause  any  numerical 
problems;  however,  the  inability  of  high  water  at  the  coast  to  flow  inland  over  previously  dry  grid 
cells  can  reduce  the  accuracy  of  the  prediction  of  SSH  near  the  coast  and  does  not  allow  prediction 
of  the  extent  of  inundation  and  flooding  that  might  occur  (Oey  et  al.  2007). 

If  the  depth  of  the  water  at  a  shallow  grid  cell  decreases  to  the  extent  that  the  grid  cell  drys  out, 
i.e.,  the  SSH  drops  to  or  below  the  ocean  bottom  at  that  grid  cell,  this  currently  causes  a  numerical 
instability  in  NCOM  that  will  terminate  the  run.  Hence,  without  a  WAD  capability,  one  is  faced 
with  setting  a  minimum  depth  at  grid  cells  in  shallow  areas  near  the  coast  that  is  sufficiently  deep 
that  the  grid  cells  will  not  dry  out  during  a  model  run.  This  involves  either  setting  a  grid  cell  that 
is  too  shallow  to  be  a  land  point,  or  increasing  the  depth  of  the  grid  cell  sufficiently  that  it  will  not 
dry  out.  Hence,  the  fidelity  of  the  bottom  depth  in  these  areas  tends  to  be  compromised,  which 
can  reduce  the  accuracy  of  the  model  predictions. 

The  maximum  drop  in  sea  level  at  points  near  a  coast  varies  considerably  and  is  mainly  affected 
by  the  amplitude  of  the  tides  and  the  amount  of  “set  down”  of  the  SSH  caused  by  the  winds.  The 
amount  of  set  up  or  set  down  at  the  coast  depends  on  the  strength  and  direction  of  the  winds  and 
the  shape  of  the  local  coastline  and  bathymetry.  For  example,  the  gradient  of  the  SSH  between  the 
mouth  and  head  of  a  bay  is  approximately  proportional  to  the  distance  between  the  mouth  and  the 
head  of  the  bay  and  the  magnitude  of  the  wind  stress  and  is  inversely  proportional  to  the  water 
depth.  Thus,  the  total  set  up  or  set  down  at  the  head  of  the  bay  increases  when  the  wind  stress 
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and  the  length  of  the  bay  increase  and  when  the  water  depth  is  reduced.  Hence,  a  strong  wind 
blowing  away  from  the  head  of  a  long,  shallow  bay  towards  the  mouth  of  the  bay  is  a  situation  that 
can  lead  to  a  significant  drop  in  the  SSH  at  the  head  of  the  bay  and  may  lead  to  drying  of  grid 
cells  near  the  head  of  the  bay  if  the  grid  cells  are  sufficiently  shallow. 

In  this  report,  a  distinction  is  made  between  “extensive”  and  “occasional”  WAD.  Extensive 
WAD  occurs  when  there  are  extensive  areas  where  WAD  is  regularly  occurring  during  the  model 
run.  This  typically  occurs  when  relatively  high  horizontal  grid  resolution  is  being  used  and  there  are 
relatively  low-lying  land  areas  and/or  shallow  sea  areas  near  the  coast  and  the  tides  are  relatively 
large  so  that  WAD  is  frequently  occurring.  A  prime  example  of  such  an  area  is  Cook  Inlet,  Alaska, 
which  has  a  tidal  range  of  over  10  m  near  its  head  and  contains  large  areas  of  tidal  mud  flats  that 
are  uncovered  during  low  tides  and  covered  up  during  high  tides  (Oey  et  al.  2007). 

Occasional  WAD  occurs  when  a  domain  does  not  have  extensive  WAD  areas  relative  to  the 
grid  resolution  being  used,  and  drying  of  grid  cells  occurs  only  occasionally,  e.g.,  due  to  the  tides 
or  to  set  down  of  the  SSH  at  the  coast  or  at  the  head  of  a  long,  shallow  bay  due  to  offshore  winds 
or  a  combination  of  both.  In  such  cases,  a  WAD  capability  allows  the  ocean  model  to  continue 
running  smoothly  during  periods  when  one  or  more  relatively  shallow  grid  cells  have  dried  out. 

This  report  is  organized  as  follows.  The  implementation  of  WAD  in  NCOM  is  discussed  in 
Section  2.  Tests  of  the  WAD  scheme  are  presented  in  Section  3.  Some  limitations  of  the  WAD 
scheme  implemented  in  NCOM  are  discussed  in  Section  4.  A  summary  of  the  report  is  presented 
in  Section  5. 

2.  IMPLEMENTATION  OF  WAD  IN  NCOM 


2.1  Scheme  used  to  implement  WAD  in  POM  by  Leo  Oey 


WAD  was  fairly  recently  implemented  in  the  Princeton  Ocean  Model  (POM)  (Blumberg  and 
Mellor  1987)  by  Leo  Oey  (Oey  2005,  2006)  using  a  relatively  simple  scheme,  and  the  WAD  imple¬ 
mented  in  NCOM  is,  in  part,  based  on  the  methods  used  by  Oey.  The  key  part  of  the  implemen¬ 
tation  of  WAD  in  POM  occurs  in  the  solution  of  POM’s  barotropic  or  free-surface  mode,  which 
involves  the  solution  of  the  equations  for  the  update  of  the  SSH  (  and  the  barotropic  transports  U 
and  V  (the  barotropic  transport  is  the  water  depth  times  the  depth-averaged  velocity).  These  are 
similar  to  the  barotropic  equations  for  a  single-layer,  shallow-water  model.  The  depth-integrated 
continuity  equation  can  be  written  as 


d(  du  dv 

dt  Ox  dy  + 


(1) 


and  the  depth-integrated  momentum  equations  can  be  written  as 


du 

~dt 


- gD ^  +  Fu, 
ox 


(2) 


dV 

~dt 


- gD ^  +  fv, 

dy 


(3) 


where  t  is  the  time,  x  and  y  are  the  horizontal  coordinates,  Q  is  a  depth-integrated  volume  source 
term,  g  is  the  acceleration  of  gravity,  D  is  the  dynamic  water  depth  (i.e.,  which  changes  in  time 
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due  to  changes  in  the  SSH),  and  Fu  and  Fv  represent  the  vertical  integral  of  the  forcing  terms  from 
the  baroclinic  momentum  equations  being  applied  to  the  barotropic  equations. 

At  this  point  it  is  useful  to  define  the  static  bottom  depth  H,  which  is  the  depth  from  the  static 
(fixed)  surface  at  z  =  0  to  the  sea  bottom.  This  static  bottom  depth  H  is  defined  to  be  positive 
upward  in  both  POM  and  NCOM,  hence,  its  value  is  negative  at  sea  points  where  the  bottom  lies 
below  z  =  0  and  is  positive  at  points  where  the  bottom  lies  above  z  =  0.  The  dynamic  water 
depth  D  is  the  total  thickness  of  the  water  column  between  the  bottom  and  the  free  surface  </.  The 
dynamic  water  depth  is  computed  as  D  =  £  —  H  and  should  always  be  >  0. 

POM  solves  its  barotropic  equations  using  what  is  referred  to  as  a  split-explicit  scheme  whereby 
these  equations  are  solved  explicitly  with  a  smaller  time  step  than  the  rest  of  POM’s  equations,  i.e. , 
the  equations  for  the  baroclinic  (depth-dependent)  velocities  and  the  temperature  and  salinity  fields. 
Note  that  POM  also  includes  barotropic  (i.e.,  depth-integrated)  forms  of  the  horizontal  advection 
and  Coriolis  terms  in  its  depth-integrated  momentum  equations,  which  are  updated  during  the  time 
integration  of  the  barotropic  mode.  Since  these  terms  are  not  essential  for  the  discussion  presented 
here,  they  have  not  been  explicitly  included  in  Eq.  (2)-(3). 

The  temporal  numerical  scheme  used  to  update  the  barotropic  equations  explicitly  in  POM, 
which  is  the  leapfrog  time-differencing  scheme,  can  be  written  for  the  depth-integrated  continuity 
equation  as 

Cn+1  =  c~l  +  2At[-(5x(AyuUn)/Ax  +  5y{AxvVn)  /  Ay)  +  Qn }  (4) 

and  for  the  depth-integrated  momentum  equations  as 

Un+ 1  =  Un~l  +  2At[-gDu6xC/Axu  +  F£]  (5) 

Vn+1  =  yn_i  +  2At[-gDvSyC/Ayv  +  F"]  (6) 

where  the  superscripts  n  +  1,  n,  and  n  —  1  denote  values  of  the  prognostic  variables  at  the  new,  cur¬ 
rent,  and  previous  time  levels  of  the  leapfrog  time-differencing  scheme,  respectively,  the  subscripts 
u  and  v  denote  values  at  the  u  and  v  points,  respectively  (both  POM  and  NCOM  use  an  Arakawa 
C-type  grid  where  the  velocities  are  defined  at  the  center  of  the  faces  of  the  main  grid  cells),  At 
denotes  the  time  step,  Ax  and  Ay  denote  grid  spacings  in  the  x  and  y  directions,  respectively,  and 
the  symbols  5X  and  Sy  denote  taking  differences  in  the  x  and  y  directions,  respectively. 

What  makes  the  update  of  these  equations  explicit  is  that  everything  that  is  needed  to  compute 
the  new  values  of  £,  U,  and  V,  which  is  on  the  right-hand-side  (RHS)  of  Eq.  (4)— (6),  is  known,  and 
so  the  new  values  can  be  easily  calculated. 

These  explicit  equations  must  be  solved  with  a  small  time  step  because,  for  explicit  equations, 
the  time  step  must  be  smaller  than  the  time  required  for  the  fastest  propagating  waves  or  signals 
described  by  the  equations  to  travel  the  length  of  a  single  grid  cell.  For  these  barotropic  equations, 
the  fastest  signals  are  the  surface  gravity  waves,  which  have  a  speed  of  csgw  =  (gD)1/2 .  In  the  deep 
ocean,  i.e.,  where  the  depth  is  on  the  order  of  4000  m  or  more,  the  speed  of  the  surface  gravity 
waves  is  greater  than  200  m/s.  By  contrast,  the  fastest  signals  in  the  baroclinic  equations  are  the 
speeds  of  the  internal  gravity  waves  and  ocean  currents,  which  are  generally  not  much  more  than 
3-4  m/s  and  are  usually  somewhat  less. 
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Oey’s  WAD  method  is  the  following:  (1)  All  grid  cells  that  are  “wettable”,  i.e.,  either  can  be 
wet  or  are  wet,  always  remain  in  the  “computational  domain”  and  are  (usually)  always  covered  by 
a  thin  film  of  water.  (Grid  cells  not  defined  to  be  in  the  computational  domain  are  considered  to  be 
land  that  can  never  be  wet.)  (2)  After  the  new  barotropic  variables  (£,  U ,  and  V)  are  computed, 
if  the  new  water  depth  D  in  any  grid  cell  is  less  than  some  minimum  specified  depth  Dmin  (e.g., 
Dmin  =  0.1  nr),  then  the  transports  directed  “out”  of  those  grid  cells  are  set  to  zero.  (3)  On  the 
next  time  step,  step  (2)  is  repeated. 


Because  the  new  transports  computed  by  Eq.  (5)-(6)  will  be  used  to  update  the  SSH  on  the 
next  time  step  using  Eq.  (4),  Oey’s  step  (2)  insures  that  at  grid  cells  with  a  water  depth  less  than 
Dmin,  the  water  depth  will  not  drop  any  further. 


Oey’s  step  (2)  presumes  that  the  water  depth  D  in  a  grid  cell  will  not  drop  from  above  Dmin 
to  zero  or  below  in  a  single  time  step.  Because  of  the  small  time  step  used  for  the  barotropic 
equations,  and  depending  on  the  value  of  Dmin ,  this  is  usually  the  case.  However,  even  with  the 
small  time  steps  typically  used  with  POM’s  split-explicit  free  surface,  it  could  possibly  happen  that 
D  becomes  <  0.  If  this  does  happen,  then  either  Dmm  must  be  increased,  or  the  time  step  being 
used  for  the  barotropic  equations  must  be  reduced. 


A  situation  that  increases  the  likelihood  of  the  water  depth  in  a  grid  cell  dropping  more  than 
Dmin  and  drying  out  in  a  single  barotropic  time  step  is  having  a  shallow  grid  cell  that  can  dry  out 
next  to  a  relatively  deep  grid  cell.  On  a  typical  sigma  coordinate  grid,  the  bottom  depth  at  the 
velocity  point  on  the  interface  between  two  grid  cells  is  the  mean  of  the  bottom  depth  of  the  two 
grid  cells.  Hence,  if  one  of  the  grid  cells  is  relatively  deep,  the  depth  at  the  velocity  point  will  also 
be  relatively  deep,  and  the  transport  between  the  two  grid  cells,  which  is  proportional  to  the  area 
of  the  grid-cell  face  that  the  two  grid  cells  share,  could  potentially  be  fairly  large  and,  hence,  more 
likely  to  dry  out  the  shallow  grid  cell  during  a  time  step. 


One  strategy  to  reduce  the  likelihood  of  this  happening  is  to  avoid  having  a  very  deep  grid  cell 
next  to  a  shallow  grid  cell  that  can  dry  out.  Another  strategy  is  to  define  the  depth  of  the  velocity 
point  between  the  two  grid  cells  to  be  equal  to  that  of  the  shallower  of  the  two  grid  cells  rather 
than  the  mean  depth  of  the  two  grid  cells.  Using  either  of  these  strategies  will  tend  to  reduce  the 
rate  of  transport  of  water  out  of  the  shallow  grid  cell  during  a  time  step  and,  hence,  tend  to  reduce 
the  rate  at  which  it  dries  out.  Note  that  if  such  a  modification  of  the  bottom  depth  at  velocity 
points  was  being  used,  it  would  only  need  to  be  done  in  shallow  areas  where  WAD  can  occur. 


An  advantage  of  the  WAD  scheme  implemented  by  Oey  (2005,  2006)  is  that  not  many  additional 
changes  to  POM  were  required  to  adapt  the  rest  of  the  model  for  WAD.  The  main  additional 
changes  that  are  required  for  WAD  are  as  follows.  (1)  The  baroclinic  velocities  are  set  to  zero  at 
the  boundaries  of  grid  cells  that  have  “dried  out”,  i.e.,  where  D  <  Dmin-  (2)  Static  bottom  depths 
H  that  are  greater  than  zero,  i.e.,  at  what  would  normally  be  land  areas,  need  to  be  allowed  for. 
(3)  Calculations  that  depend  on  the  thickness  or  depth  of  a  grid  cell,  such  as  the  bottom  drag 
coefficient  and  the  solar  extinction,  may  need  to  be  updated  more  frequently,  since  the  values  at 
shallow  grid  cells,  especially  those  that  are  beginning  to  flood  or  are  near  to  drying,  can  undergo 
relatively  large  changes  in  a  short  period  of  time. 
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2.2  Implementation  of  WAD  in  NCOM’s  Explicit  Barotropic  Solver 

The  WAD  scheme  of  Oey  (2005,  2006),  described  in  the  previous  section,  was  implemented  in 
NCOM’s  explicit  barotropic  solver  for  the  new  SSH  and  barotropic  transports.  This  routine  is  used 
to  provide  a  “preliminary”  estimate  of  the  new  SSH  when  NCOM  is  run  in  the  usual  way  using 
an  implicit  solution  for  the  new  SSH,  and  is  also  used  to  compute  the  new  SSH  and  barotropic 
transports  explicitly  when  NCOM  is  run  in  fully  explicit  mode  using  the  same  (small)  time  step 
for  both  the  barotropic  and  baroclinic  equations. 

Note  that  NCOM  is  not  usually  run  in  fully  explicit  mode  because,  using  the  small  time  step 
that  is  needed  to  explicitly  solve  the  barotropic  equations  to  solve  the  baroclinic  equations  results 
in  NCOM  taking  a  great  deal  of  time  to  run,  i.e. ,  the  size  of  the  time  step  required  to  solve  the 
barotropic  equations  explicitly  can  be  more  than  50  times  smaller  than  that  required  to  solve  the 
baroclinic  equations  if  the  model  domain  contains  deep  areas.  However,  this  option  is  provided  for 
testing  purposes,  or  in  case  a  fully  explicit  solution  is  desired,  e.g.,  for  comparison  with  the  implicit 
solution  for  the  new  SSH. 

2.3  Implementation  of  WAD  in  NCOM’s  Implicit  Barotropic  Solver 

To  implement  WAD  with  NCOM’s  implicit  solution  of  the  barotropic  equations  (Martin  2000), 
a  variation  of  Oey’s  WAD  scheme  is  used.  After  the  new  SSH  is  computed  (using  an  iterative 
solver),  the  new  water  depth  in  each  grid  cell  is  inspected  and,  for  grid  cells  whose  depths  have 
fallen  below  a  prescribed  minimum  Dmin,  the  volume  fluxes  at  cell  faces  that  are  directed  “out”  of 
the  drying  grid  cells  are  set  to  zero  by  setting  the  solver  coefficient  for  those  grid-cell  faces  to  zero, 
and  then  the  free-surface  solver  is  rerun. 

This  iterative  procedure  is  repeated  until  convergence  occurs,  i.e.,  until  no  grid-cell  depths  drop 
below  Drnjn  or  the  maximum  difference  in  the  SSH  from  the  calculation  of  the  SSH  on  the  previous 
iteration  falls  below  a  small  prescribed  value  (currently  set  to  10~6  m).  Note  that  volume  fluxes  set 
to  zero  on  earlier  iterations  during  the  current  time  step  remain  set  to  zero  during  later  iterations. 
However,  as  with  Leo  Oey’s  scheme  in  POM,  everything  starts  over  again  on  the  next  time  step, 
i.e.,  all  the  grid  cells  in  the  “computational  domain”  are  included  in  the  initial  calculation  of  the 
new  SSH  on  the  next  time  step.  This  allows  “dry”  grid  cells  to  readily  flood  when  the  SSH  in  the 
adjoining  grid  cells  begins  to  rise. 

Note  that,  with  this  scheme,  the  water  depth  D  at  grid  cells  within  the  computational  domain 
is  not  allowed  to  fall  much  below  the  minimum  specified  depth  Dmin.  For  this  reason,  the  scheme 
has  proven  to  be  fairly  robust  in  the  testing  done  to  date. 

Since  the  water  depth  D  at  grid  cells  within  the  computational  domain  cannot  fall  much  below 
the  value  specified  for  Dmini  one  could  consider  making  Drnin  fairly  small,  so  that  the  layer  of  water 
left  behind  in  “dry”  areas  will  tend  to  be  relatively  shallow.  However,  a  drawback  of  this  is  that,  as 
the  water  depths  become  very  shallow,  larger  velocities  may  occur  in  the  areas  where  the  grid  cells 
are  shallow,  which  would  increase  the  possibility  of  violating  the  Courant-Friedrichs-Lewy  (CFL) 
limitation  for  explicit  numerical  advection,  which  is  that  the  advection  of  a  field  cannot  exceed  the 
width  (or  thickness  for  vertical  advection)  of  a  grid  cell  in  a  single  time  step.  We  typically  set  the 
value  for  Dmin  in  our  coastal  simulations  to  10  cm.  Even  with  this  value  of  Dmin,  the  time  step  for 
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simulations  that  include  WAD  sometimes  has  to  be  reduced  to  prevent  violation  of  the  advective 
CFL  limitation  in  the  WAD  areas. 

Note  that  reduction  of  the  time  step  is  not  always  needed  when  WAD  is  added,  or  when  using 
WAD  and  the  value  of  Drnin  is  reduced,  but  if  the  time  step  does  have  to  be  reduced,  it  can  result 
in  significant  additional  cost  for  the  simulation,  e.g.,  if  the  time  step  is  halved,  the  time  required 
for  the  simulation  will  be  almost  doubled.  This  is  significantly  more  than  the  cost  of  just  the 
additional  calculations  required  for  WAD,  which  is  typically  less  than  20%.  We  have  conducted 
WAD  simulations  with  values  of  Dmin  of  10,  5,  2,  1,  and  0.1  cm,  and  these  have  run  robustly  as  long 
as  the  time  step  was  reduced  sufficiently  to  avoid  exceeding  the  advective  CFL  limitation.  Finding 
the  “maximum”  time  step  that  will  work  for  a  given  domain  or  a  particular  situation  usually  takes 
some  trial  and  error,  though  this  is  also  the  case  for  simulations  without  WAD. 

2.4  General  Modifications  of  NCOM  for  WAD 

At  grid-cell  faces  at  which  the  barotropic  transport  is  set  to  zero  during  the  calculation  of  the 
barotropic  equations  (i.e.,  to  prevent  the  water  depth  D  at  a  grid  cell  from  dropping  below  Dmin ), 
the  baroclinic  velocities  are  also  set  to  zero. 

Without  WAD,  the  static  bottom  depth  H  at  active  (sea)  grid  cells  in  the  computational 
domain  is  always  less  than  zero.  However,  with  WAD  the  static  bottom  depth  at  an  active  grid 
cell  can  be  set  to  be  >  0.  Hence,  locations  in  the  NCOM  computer  code  where  the  use  of  a  value 
of  H  >  0  would  generate  unphysical  values  must  be  modified.  This  is  typically  done  by  setting  a 
maximum  allowable  value  for  H  that  is  less  than  zero  (e.g.,  -0.1  m)  to  be  used  in  such  calculations. 

Some  parameterizations  in  NCOM  that  depend  on  the  layer  depth  or  thickness,  such  as  the 
calculation  of  the  the  solar  extinction  and  the  bottom  drag  coefficient,  are  usually  computed  only 
at  the  start  of  a  model  run  based  on  the  static  layer  depths  or  thicknesses.  However,  with  WAD, 
the  calculation  of  the  solar  extinction  and  the  bottom  drag  coefficient  with  the  static  layer  depths 
or  thicknesses  will  be  inaccurate  and/or  ill-defined  in  WAD  areas,  and  both  the  solar  extinction  and 
the  bottom  drag  coefficient  can  change  significantly  at  grid  cells  undergoing  WAD.  Hence,  options 
are  provided  to  update  these  quantities  each  time  step  using  the  current  dynamic  layer  depths 
and  thicknesses,  which  change  on  the  sigma-coordinate  grid  with  changes  in  the  SSH.  Note  that, 
previously,  calculation  of  the  bottom  drag  coefficient  each  time  step  using  the  dynamic  bottom-layer 
thickness  was  provided  when  ocean-wave  coupling  was  being  used  to  allow  for  the  enhancement  of 
the  bottom  drag  coefficient  by  the  wave  motion  near  the  bottom  (i.e.,  in  the  wave-bottom  boundary 
layer),  which  has  a  large  variation  in  time  as  the  waves  change. 

For  WAD  a  new  input  flag,  indwad ,  is  provided  for  NCOM  to  define  whether  WAD  will  be 
computed  during  the  run.  This  flag  is  set  to  zero  if  there  is  no  WAD  (this  is  the  default  value)  and 
must  be  set  to  a  value  of  one  if  WAD  is  to  be  accounted  for.  Another  new  input  flag,  iwadiag,  can 
be  used  to  increase  the  amount  of  diagnostic  printout  associated  with  the  WAD.  The  default  value 
of  iwadiag  is  zero,  and  larger  values  increase  the  amount  of  diagnostic  output.  This  output  can  be 
used  to  help  diagnose  any  problems  associated  with  the  WAD  and  to  check  the  conservation  of  the 
WAD  scheme. 

The  value  of  Dmin  can  be  specified  as  an  input  to  NCOM  if  a  value  other  than  the  default 
value  of  0.1  m  is  desired.  However,  note  that  within  the  NCOM  code  and  within  NCOM’s  input 
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parameter  files,  this  parameter  is  called  Dwet,  which  was  considered  to  be  a  more  unique  and 
representative  name. 

When  computing  the  new  value  of  the  SSH  and  checking  for  values  of  D  <  Dmin  for  WAD,  it 
is  not  necessary  to  check  points  where  the  value  of  D  is  relatively  large  that  are  unlikely  to  dry 
out  during  that  time  step.  Hence,  the  first  time  this  check  is  being  conducted  on  a  particular  time 
step,  NCOM  makes  a  list  of  the  grid  cells  within  the  computational  domain  where  D  is  less  than 
a  specified  minimum  value  DV)acp  and  only  these  points  are  checked  for  D  <  Dmin  on  that  time 
step.  This  reduces  the  amount  of  time  that  NCOM  spends  doing  this  checking.  The  default  value 
of  Dwad  is  set  to  4  nr,  but  a  different  value  can  be  specified  as  an  input  parameter.  Experience  has 
shown  that  Dwad  needs  to  be  set  to  a  value  on  the  order  of  the  maximum  tidal  amplitude  within 
the  domain,  or  on  the  order  of  the  maximum  wind  setup  or  set  down,  whichever  is  greater.  For 
example,  a  value  of  Dwaci  =  3  nr  was  found  to  be  adequate  for  simulations  of  the  tide  in  Cook 
Inlet  where  the  maximum  tidal  amplitude  is  about  5  nr.  If  Dwad  is  set  too  small,  the  calculation 
of  the  new  SSH  may  not  converge  and  the  NCOM  simulation  will  terminate.  Hence,  if  an  NCOM 
simulation  terminates  during  the  calculation  of  WAD  during  the  update  of  the  SSH,  it  may  be  that 
the  value  of  Dwad  being  used  is  too  small. 

When  WAD  is  being  used,  the  surface  values  of  the  land-sea  mask,  which  define  which  grid  cells 
are  inside  and  which  are  outside  the  computational  domain,  can  no  longer  easily  be  determined 
from  the  two-dinrensional  (2D)  array  of  static  bottom  depths  H.  For  this  reason,  a  2D,  surface, 
land-sea  mask  is  now  added  to  the  NCOM  input  horizontal  grid  file  during  the  setup  of  an  NCOM 
run.  This  2D  land-sea  mask  becomes  the  seventh  2D  record  within  the  horizontal  grid  Hie,  and 
follows  the  six  previous  2D  arrays  in  this  file  containing:  (1)  longitude,  (2)  latitude,  (3)  grid  spacing 
in  x,  (4)  grid  spacing  in  y,  (5)  static  bottom  depth  H  at  sea  (i.e.,  all  wettable)  points,  and  (6)  the 
angle  of  the  grid  with  respect  to  the  local  longitude  and  latitude.  When  NCOM  starts  up,  the 
land-sea  mask  is  read  from  the  horizontal  grid  file  if  WAD  is  being  used.  If  WAD  is  not  being  used, 
the  land-sea  mask  is  computed  from  the  bottom  depth  array  as  was  done  before  without  WAD,  i.e., 
points  with  H  <  0  are  defined  as  sea  points  and  are  inside  the  computational  domain  and  all  others 
are  defined  as  land  points  and  are  outside  the  computational  domain.  This  (and  the  default  value 
of  input  parameter  indwad  being  zero)  helps  maintain  backward  compatability  for  the  updated 
NCOM  code  that  includes  WAD,  with  older  sets  of  NCOM  input  files  that  were  set  up  before  the 
WAD  was  implemented. 

2.5  Bottom  drag 

The  bottom  drag  plays  an  important  role  in  the  dynamics  of  WAD,  i.e.,  very  shallow  flows  in 
WAD  areas  are  often  mainly  determined  by  a  balance  between  the  surface  pressure  gradient  due 
to  the  slope  of  the  SSH  and  the  bottom  drag.  Hence,  the  calculation  of  the  bottom  drag  and  the 
specification  of  the  bottom  drag  coefficient  is  an  important  aspect  of  WAD  simulations.  For  these 
reasons,  the  bottom  drag  parameterization  used  in  NCOM  will  be  discussed  here. 

NCOM  uses  a  quadratic  parameterization  of  the  bottom  drag  (Martin  2000);  hence,  the  bottom 
stress  terms  on  the  RHS  of  NCOM’s  u  and  v  momentum  equations  are  —  Cbu\v\  and  —Cyv |v|, 
respectively,  where  65  is  the  bottom  drag  coefficient,  and  u  and  v  are  the  two  components  of  and 
|  v |  is  the  magnitude  of  the  horizontal  velocity  in  the  bottom  layer. 

The  value  of  the  bottom  drag  coefficient  Cb  used  in  NCOM  can  be  specified,  or  can  be  calculated 
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in  terms  of  the  bottom  layer  thickness  A zb  and  the  bottom  roughness  zQ  as 
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where  k  =  0.4  is  von  Karman’s  constant,  and  Cb  .  is  a  minimum  allowed  value  for  Cb-  If  the 
bottom  roughness  z0  is  set  to  zero,  the  bottom  drag  is  just  set  equal  to  Cbmin.  This  expression 
for  Cb  assumes  (i.e.,  is  derived  assuming)  a  logarithmic  boundary  layer  velocity  profile  near  the 
bottom.  For  a  flow  being  treated  as  a  single  layer,  the  value  of  A zb  in  Eq.  (7)  is  just  the  water 
depth. 


A  commonly  used  alternative  for  computing  the  drag  on  a  flow,  which  is  frequently  used 
in  engineering,  is  the  Manning  formulation.  With  the  Manning  formulation,  the  bottom  drag 
coefficient  for  a  shallow-water  flow  that  is  being  treated  as  a  single-layer  flow  is  computed  as 

Cb  =  gn2/D V3  (8) 

where  g  is  the  gravitational  constant,  D  is  the  water  depth,  and  n  is  the  Manning  roughness 
coefficient,  which  depends  on  the  characteristics,  e.g.,  the  roughness,  of  the  bottom.  Since  the 
bottom  drag  coefficient  Cb  is  dimensionless,  the  Manning  roughness  coefficient  has  units  of  s/m1/3. 
Values  of  the  Manning  roughness  coefficient  can  be  found  in  the  literature  for  a  wide  range  of 
material  and  bottom  types. 


The  use  of  the  Manning  formulation  for  the  bottom  drag  is  not  currently  provided  as  an  option 
in  NCOM.  However,  in  both  of  the  bottom  drag  coefficient  formulations  described  by  Eq.  (7) 
and  (8),  the  bottom  drag  coefficient  increases  as  the  water  depth  decreases,  though  at  somewhat 
difference  rates.  Hence,  when  describing  a  single-layer  flow,  an  approximate  correspondence  can 
be  found  between  the  bottom  roughness  za  and  the  Manning  friction  coefficient  n  for  a  specific, 
limited  range  of  water  depths. 


2.6  Pre-  and  Post-Processing  Modifications  for  WAD 

Calculations  used  in  pre-processing  (i.e.,  setting  up)  an  NCOM  simulation  as  well  as  post¬ 
processing  of  NCOM  output  must  allow  for  the  use  of  WAD  if  there  are  values  of  the  static  bottom 
depth  H  >  0  within  the  computational  domain. 

Setting  up  the  bathymetry  for  a  domain  involves  a  number  of  procedures  that  are  sensitive  to 
the  value  of  H  being  >  0.  These  procedures  include:  calculation  of  the  land-sea  mask,  smooth¬ 
ing/filtering  of  the  bathymetry,  calculation  of  steep  bottom  slopes,  reduction  of  steep  bottom  slopes, 
and  the  determination  of  a  single,  main,  contiguous,  computational  domain. 

For  all  of  these  procedures,  values  of  H  >  0  can  be  accommodated  by  defining  (a)  a  minimum 
elevation  for  land  points  that  are  defined  to  be  outside  the  computational  domain  and/or  (b)  a 
reference  height  that  is  above  the  highest  allowable  elevation  of  grid  cells  that  are  defined  to  be 
within  the  computational  domain.  Setting  a  minimum  elevation  to  define  land  points  that  are 
outside  the  computational  domain  provides  a  simple  rule  for  computing  the  land-sea  mask  that 
defines  the  computational  domain.  Alternatively,  the  rule(s)  used  for  defining  the  computational 
domain  can  vary  spatially,  but  this  is  more  complex  to  implement.  Defining  a  reference  height 
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for  the  grid  cells  within  the  computational  domain  is  useful  for  providing  a  reference  height  for 
computing  and  adjusting  bottom  slopes.  Note  that  smoothing  bottom  depths  and  reducing  steep 
slopes  within  and  adjacent  to  WAD  areas  can  increase  the  robustness  of  the  WAD  calculations 
within  an  NCOM  simulation  and  may  reduce  the  need  to  decrease  the  time  step  to  avoid  violating 
the  advective  CFL  constraints. 

As  noted  previously,  when  providing  for  WAD,  NCOM  tries  to  maintain  a  minimum  water 
depth  Dmin  at  all  grid  cells  within  the  computational  domain.  Consequently,  when  the  SSH  (£)  is 
initialized  in  an  NCOM  setup  program,  the  SSH  at  grid  cells  within  the  computational  domain  can 
be  defined  so  that  the  initial  dynamic  water  depth  D  =  £  —  H  is  >  Dmin.  However,  NCOM  can 
adjust  if  a  different  value  of  Dmin  was  used  to  define  the  initial  SSH  as  long  as  D  >  0. 

For  convenience,  pre-  and  post-processing  of  NCOM  input  and  output  frequently  involves 
calculating  a  three-dimensional  (3D)  array  of  static  or  dynamic  grid-cell  interface  depths.  Such  an 
array  can  be  used  to  define  the  vertical  structure  of  most  grids,  including  NCOM’s  sigma,  sigma/z- 
level,  and  generalized  vertical  coordinate  grids,  and  is  useful  for  computing  or  accessing  depth- 
dependent  quantities  and  values  from  data  bases  for  model  setup  and  for  processing  model  output. 
In  the  case  where  WAD  is  being  used  and  there  are  values  of  H  >0,  it  is  usually  sufficiently  accurate 
for  most  pre-  and  post-processing  needs  to  define  a  3D  array  of  static  interface  depths  in  which  the 
maximum  (shallowest)  layer  depths  are  defined  to  be  small,  negative  values  increasing  (however 
slightly)  upwards  towards  a  static  free  surface  at  zero  elevation,  so  that  the  depths  computed  from 
this  array  of  interface  depths  will  always  be  negative  and  the  computed  layer  thicknesses  will  always 
be  positive  (i.e.,  as  for  the  case  when  WAD  is  not  being  used).  This  is  because  many  of  the  data 
bases  and  procedures  used  in  pre-  and  post-processing  are  expecting  a  value  of  the  depth  at  a  sea 
point  that  is  <  0. 

3.  WAD  TEST  CASES 

3.1  Initial  Gaussian  Bump  Symmetry  Test 

NCOM  is  written  so  as  to  be  able  to  maintain  perfect  symmetry  for  horizontally  symmetric 
problems.  This  is  accomplished  by  grouping  terms  in  the  calculations  in  NCOM  to  maintain 
horizontal  symmetry,  e.g.,  if  differencing  the  horizontal  fluxes  on  the  six  faces  of  a  grid  cell,  the 
separate  differences  in  x,  y,  and  z  are  computed  first,  then  the  x  and  y  flux  differences  are  combined, 
and  then  the  flux  difference  in  z  is  added.  Such  grouping  of  the  calculations  makes  use  of  the  fact 
that,  in  most  computers,  the  sum  or  product  of  two  quantities  gives  the  same  result  independent 
of  which  quantity  appears  first  in  the  compute  instruction.  However,  the  sum  or  product  of  three 
or  more  quantities  will  vary,  due  to  roundoff  error,  depending  on  the  order  in  which  the  calculation 
is  performed. 

Checks  for  the  horizontal  symmetry  of  the  NCOM  solutions  for  horizontally  symmetric  prob¬ 
lems  provides  a  very  useful  check  on  the  coding  and  is  especially  useful  for  detecting  loop-indexing 
errors.  This  is  especially  useful  for  the  NCOM  code,  since  the  bulk  of  the  calculations  in  NCOM 
are  done  in  a  horizontally  asymmetric  fashion,  i.e.,  the  update  of  the  3D  momentum  fields  and 
the  3D  scalar  and  turbulence  fields  is  done  by  proceeding  through  the  model  domain  in  single  x-z 
slices.  This  is  done  to  improve  the  efficiency  of  the  use  of  high-speed  cache  memory  when  there 
are  a  lot  of  grid  points  being  computed  on  each  processor.  By  only  computing  a  single,  2D,  x-z 
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slice  at  a  time  (rather  than  the  full  3D  grid),  there  is  greater  likelihood  that  more  calculations  can 
be  done  with  the  variables  being  held  within  the  high-speed  cache  before  they  are  flushed  out  by 
the  need  to  access  variables  not  currently  in  the  high-speed  cache.  In  summary,  the  asymmetric 
way  in  which  most  of  NCOM’s  updates  proceed  through  the  model  domain  makes  for  efficient  use 
of  high-speed  cache,  but  provides  a  lot  of  potential  for  coding  errors.  Hence,  a  test  for  horizontal 
symmetry  is  usually  the  first  test  that  is  performed  on  NCOM  after  extensive  coding  changes  have 
been  made. 

The  typical  horizontal  symmetry  problem  that  is  used  is  set  up  within  a  square  domain.  The 
lateral  boundaries  may  be  open  or  closed  (usually,  both  of  these  situations  are  tested).  The  initial 
condition  is  taken  to  be  a  Gaussian-shaped  bump  in  the  SSH,  centered  within  the  domain.  A 
similar  type  of  symmetric  bump  is  usually  incorporated  in  the  initial  salinity  held.  Hence,  when 
the  simulation  is  started,  both  the  SSH  bump  and  the  internal  salinity  bump  propagate  outwards 
towards  the  edges  of  the  domain  as  external  and  internal  gravity  waves,  and  either  pass  through  the 
open  boundaries,  or  are  reflected  at  the  boundaries  if  they  are  closed.  Typically,  as  much  variability 
as  possible  is  added  to  this  problem  in  order  to  make  the  test  as  thorough  as  possible,  the  only 
limitation  being  that  what  is  added  is  added  in  a  horizontally  symmetric  way,  e.g.,  spatially- variable 
horizontal  grid  spacing,  spatially-variable  bathymetry,  atmospheric  forcing,  and  river  inflows.  The 
resulting  solutions  should  have  an  eight-fold  symmetry  when  the  Coriolis  term  is  taken  to  be  zero, 
and  a  four-fold  symmetry  if  the  Coriolis  term  is  taken  to  be  constant.  NCOM  provides  subroutines 
that  automatically  check  for  these  symmetries  when  this  test  is  being  conducted. 

To  test  NCOM’s  WAD  scheme,  WAD  areas  were  added  to  the  grid  in  a  symmetric  fashion  to 
be  exposed  or  submerged  as  the  initial  Gaussian  SSH  bump  propagates  outward  from  the  center 
of  the  grid.  Note  that  our  initial  symmetry  tests  of  the  WAD  scheme  did  expose  some  problems 
with  the  changes  that  had  been  made  to  incorporate  the  WAD  (i.e.,  as  it  has  many  times  in  the 
past,  the  symmetry  test  proved  its  usefulness).  After  correction  of  these  problems,  the  solutions 
were  symmetric. 

Note  that  in  recent  years,  we  have  found  that  successful  conduct  of  the  symmetry  tests  requires 
that  the  NCOM  code  be  compiled  at  a  fairly  low  level  of  optimization.  It  appears  that,  at  high  opti¬ 
mization  levels,  the  parentheses  employed  in  the  NCOM  code  to  specify  the  order  of  the  arithmetic 
operations  are  sometimes  ignored/over-ruled  by  the  aggressive  optimization  of  the  compiler.  When 
this  happens,  the  solutions  still  tend  to  be  “correct”,  i.e.,  within  the  limits  of  computer  roundoff 
error,  but  may  not  be  perfectly  symmetric. 

3.2  Sloshing  around  a  Parabolic-Shaped  Basin 

This  test  problem  was  used  by  Guan  et  al.  (2013)  to  test  a  WAD  scheme.  Note  that  this  test 
problem  was  originally  described  by  Thacker  (1981).  The  test  involves  the  propagation  of  a  planar 
surface  wave  around  in  a  parabolic-shaped  basin.  WAD  occurs  at  the  edge  of  the  basin  as  the 
surface  wave  propagates  along  the  edge  of  the  basin.  The  shape  of  the  basin  is  given  by 

h(x,y)  =  ^(x2  +  y2),  (9) 

az 

where  h  is  the  bottom  depth  relative  to  a  zero  reference  level,  ho  =  — 10  m  is  the  water  depth  at 
the  center  of  the  basin,  a  =  8025.5  m  is  the  radius  of  the  basin  at  zero  elevation,  and  x  and  y  are 
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the  horizontal  coordinates.  Note  that  h  was  defined  to  be  positive  downward  in  Guan  et  al.  (2013), 

but  has  been  defined  to  be  positive  upward  here  to  be  consistent  with  the  sign  convention  used  in 

NCOM. 

This  problem  has  an  exact  analytical  solution  given  by 

C(x,y,t)  = - 2^[2xcos(o;t)  +  2ysin(u;f)  —  A]  —  ho,  (10) 

u(t)  =  —  Ansin(wt),  (11) 

v(t)  =  — Aiu  cos(aA),  (12) 

where  (  is  the  SSH,  A  =  a/10  =  802.55  m  is  a  constant  that  determines  the  amplitude  of  the 
motion,  u  =  (2gho)0'5  /  a  =  2ir/T  is  the  rotational  frequency,  g  =  9.81  m2/s  is  the  acceleration  of 
gravity,  T  =  3600  s  is  the  period  of  rotation,  t  is  the  time  in  s,  and  u  and  v  are  the  components  of 
the  velocity  in  the  x  and  y  directions,  respectively.  Note  that  the  velocity  is  always  at  a  right  angle 
to  the  slope  of  the  SSH,  and  the  sign  used  in  Eq.  (12)  for  the  v  velocity  determines  the  direction 
of  rotation,  which  will  be  clockwise  in  this  case. 

As  in  Guan  et  al.  (2013),  the  problem  is  solved  in  a  20  x  20  km  domain  with  a  200  x  200  grid 
and  a  uniform  horizontal  grid  resolution  of  100  m.  Only  one  layer  is  used  in  the  vertical.  The  initial 
SSH  and  velocities  are  as  described  by  Eq.  (10)-(12)  at  time  t  =  0.  The  bottom  drag  is  set  to 
zero.  The  time  step  is  10  s.  The  third-order  upwind  scheme  (Holland  1998)  is  used  for  momentum 
advection  in  the  NCOM  simulation. 

We  initially  ran  the  simulation  with  a  minimum  water  depth  Dmin  within  the  computational 
domain  of  0.1  m.  However,  the  rotation  speed,  which  should  be  one  revolution  per  hour,  was  about 
0.8%  fast.  Reducing  the  minimum  water  depth  Drnin  to  0.01  m  reduced  the  rotation  speed  error  to 
about  0.07%.  We  also  initially  ran  with  a  small  but  non  zero  value  of  the  bottom  drag  coefficient, 
with  a  maximum  value  of  0.000042;  however,  this  resulted  in  damping  of  the  solution  of  about 
2%  per  rotation.  Setting  the  bottom  drag  to  zero,  which  it  should  be  for  this  idealized  problem, 
noticeably  reduced  the  damping.  No  reduction  of  the  time  step  was  needed  to  incorporate  these 
changes,  though  reducing  the  value  of  Dmin  from  0.1  to  0.01  m  increased  the  maximum  value  of  the 
horizontal  advective  CFL  values  from  0.31  to  0.45  (the  numerical  advection  is  stable  for  maximum 
advective  CFL  values  below  one). 

Figure  1  shows  a  comparison  of  the  SSH  simulated  by  the  model  with  the  SSH  computed  from 
the  analytical  solution  in  Eq.  (10)  along  a  line  through  the  middle  of  the  basin  at  y  =  0  at  times  of 
1,  5.2,  5.5,  and  6  h  (the  period  of  rotation  is  1  h).  The  model  solution  still  shows  good  agreement 
with  the  analytical  solution  after  6  h. 

Some  aspects  of  the  implementation  of  WAD  illustrated  by  this  test  case  are:  conservation  of 
volume,  accurate  simulation  of  the  speed  of  rotation  of  the  wave  around  the  basin,  and  negligible 
damping  of  the  amplitude  of  the  wave  (when  the  bottom  drag  is  set  to  zero). 
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Fig.  1  Comparison  of  model-simulated  SSH  (dots)  with  analytical  solution  (dashed  line)  at  y  =  0  after  (a)  1  period, 

(b)  5.2  periods,  (c)  5.5  periods,  and  (d)  6  periods. 
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gate 


Fig.  2  —  Setup  and  dimensions  for  laboratory  experiment  for  a  dam-break  flow  over  a  triangular  bump. 

3.3  Dam  Break  Flow  over  a  Triangular  Bump 

This  test  case  involves  the  simulation  of  a  WAD  laboratory  experiment  that  was  conducted 
by  the  Universite  Libre  de  Bruxelles  and  the  Laboratoire  de  Recherches  Hydrauliques  (Chatlet), 
Belgium  under  the  supervision  of  Professor  J.M.  Hiver  (Alcruco  and  Soares  Frazao  1999).  The 
experiment  was  conducted  to  provide  data  for  testing  WAD  in  dam-break  and  inundation  numerical 
models.  The  results  from  this  laboratory  experiment  have  been  used  for  testing  numerical  models 
by  a  number  of  investigators,  e.g.,  Alcrudo  (1999),  Benkhaldoun  et  al.  (1999),  Garcia-Navarro  and 
Brufau  (1999),  Brufau  et  al.  (2002),  Cozzolino  et  al.  (2006),  Loukili  and  Soulaimani  (2007),  Liang 
and  Marche  (2009),  Singh  et  al.  (2011),  Guan  et  al.  (2013),  and  Khan  and  Lai  (2014). 

The  setup  of  the  experiment  is  illustrated  in  Fig.  2.  Water  contained  within  a  holding  tank  on 
the  left  end  of  a  channel  is  suddenly  released  by  opening  a  gate  at  the  front  (i.e.,  on  the  right  end) 
of  the  tank.  After  it  is  released,  the  water  flows  down  the  channel  and  some  of  the  water  flows  over 
a  triangular-shaped  bump.  During  the  experiment,  the  surface  elevation  was  measured  at  seven 
locations  along  the  channel.  These  SSH  measurement  locations  are  shown  in  Fig.  2  (denoted  by 
the  purple  lines)  and  are  at  distances  of  2,  4,  8,  10,  11,  13,  and  20  m  from  the  front  of  the  tank. 

NCOM  was  set  up  for  this  test  case  with  a  horizontal  grid  resolution  of  0.04  m.  Beyond  the 
end  of  the  experimental  setup  shown  in  Fig.  2,  a  catch  basin  of  length  5  m  was  used  to  collect 
the  water  flowing  out  the  right  end  of  the  channel.  Hence,  the  total  length  of  the  domain  used 
for  the  NCOM  simulation,  including  the  catch  basin,  was  43  m.  Only  3  points  were  used  in  the 
cross-channel  direction,  with  just  the  middle  point  being  a  water  point  and  the  first  and  third  points 
being  land  outside  the  computational  domain  with  free-slip  land-sea  boundaries.  (Alternatively, 
periodic  boundaries  could  have  been  used  in  the  cross-channel  direction,  in  which  case  the  results 
would  have  been  the  same,  or  the  true  width  of  the  channel  in  the  y-direction  along  with  the  side 
boundaries  could  have  been  simulated.)  Only  a  single  layer  was  used  in  the  vertical.  The  third- 
order  upwind  scheme  was  used  for  momentum  advection.  Quadratic  bottom  drag  was  used  with  a 
bottom  roughness  of  3.5  x  10-'5  m  and  a  minimum  value  of  the  bottom  drag  coefficient  of  0.0001  (see 
Section  2.5).  This  bottom  roughness  corresponds  approximately  to  a  Manning  roughness  coefficient 
of  0.0125  s/m1/3,  which  is  the  value  specified  for  the  experimental  setup  by  Alcrudo  and  Soares 
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Frazao  (1999).  These  give  values  for  the  bottom  drag  coefficient  of  0.0023  for  a  water  depth  of 
0.3  m.  Runs  were  made  with  the  minimum  water  depth  Dmin  set  to  0.1,  0.01,  and  0.001  m.  Using 
smaller  values  of  reduced  the  speed  of  the  wetting  front  and  improved  agreement  with  the 

observations.  The  results  shown  are  for  Dmin  =  0.001  m.  The  time  step  used  was  0.002  s  and 
output  from  the  model  simulation  was  saved  at  0.1-s  intervals. 

Figure  3  shows  the  SSH  from  the  model  simulation  at  5-s  intervals  from  0  to  85  s.  At  5  s, 
the  released  water  has  reached  the  top  of  the  triangular  bump.  Some  of  the  water  passes  over  the 
bump  and  some  is  reflected  back  towards  the  tank.  At  30  s,  the  wave  reflected  back  towards  the 
tank  by  the  triangular  bump  has  been  reflected  by  the  wall  on  the  left  side  of  the  tank  and  has 
started  propagating  back  towards  the  bump.  This  sequence  is  repeated  several  times,  i.e.,  some 
water  passes  over  the  triangular  bump  and  a  wave  is  reflected  by  the  bump  back  towards  the  tank 
and  is  then  reflected  by  the  back  wall  of  the  tank  back  towards  the  triangular  bump.  Three  of  these 
events  occur  during  the  first  90  s  of  the  simulation,  and  a  few  more  occur  at  later  times,  until  the 
amplitude  of  the  propagating  wave  becomes  too  small  to  get  over  the  bump. 

Figure  4  shows  a  comparison  of  the  model  simulated  SSH  with  the  measured  values  at  the  seven 
SSH  measurement  locations  shown  in  Fig.  2.  The  simulated  SSH  is  a  bit  higher  than  observed  at 
gauges  3-5  and  7,  but  the  overall  agreement  of  the  SSH  is  fairly  good  and  the  timing  of  the  SSH 
peaks  in  the  model  simulation  is  also  fairly  good. 

3.4  Dam  Break  Flow  around  a  90°  Corner 

This  test  case  involves  the  simulation  of  a  WAD  laboratory  experiment  that  was  conducted  by 
the  group  of  Professor  Y.  Zech  at  the  Universite  Catholique  de  Louvain,  Belgium  (Soares  Frazao 
and  Alcrudo  1999).  This  experiment  has  been  used  for  verification  of  WAD  models  by  other 
researchers  including  Soares  Frazao  et  al.  (1999),  Viseu  et  al.  (1999),  and  Guan  et  al.  (2013),  and 
a  later,  slightly  modified  version  of  this  experiment  was  used  by  Soares  Frazao  and  Zech  (2002) 
and  Biscarini  et  al.  (2007).  The  setup  for  this  experiment  is  shown  in  Fig.  5.  A  tank  of  water  of 
horizontal  size  2.39  by  2.44  m  is  filled  to  a  depth  of  0.53  m  with  water.  The  tank  is  connected  by  a 
gate  to  a  channel  of  width  0.495  m.  At  a  distance  of  about  4  m  from  the  tank,  the  channel  makes 
a  right-angle  turn  to  the  left.  Note  that  the  floor  of  the  channel  is  0.33  m  above  the  bottom  of 
the  tank  and  is  0.20  m  below  the  initial  level  of  the  water  within  the  tank.  When  the  gate  of  the 
tank  is  opened,  the  water  flows  from  the  tank  down  the  channel  and  around  the  90°  corner  of  the 
channel.  The  water  height  within  the  tank  and  at  five  locations  along  the  channel  was  measured 
by  gauges.  The  locations  of  these  six  gauges  are  shown  in  Fig.  5.  Note  that  the  five  gauges  within 
the  channel  are  located  in  the  center  of  the  channel. 

NCOM  was  set  up  for  this  test  case  with  a  horizontal  grid  resolution  of  0.02  m.  Beyond  the 
back  end  of  the  experimental  setup  shown  in  Fig.  5,  a  catch  basin  of  1-m  width,  0.8-rn  depth,  and 
running  the  full  length  of  the  back  of  the  numerical  grid  was  used  to  collect  the  water  flowing  out 
of  the  end  of  the  channel.  The  third-order  upwind  scheme  was  used  for  momentum  advection. 
Quadratic  bottom  drag  was  used  with  a  bottom  roughness  of  6  x  10~5  nr  and  a  minimum  value  of 
the  bottom  drag  coefficient  of  0.0001  (see  Section  2.5).  This  bottom  roughness  gives  a  bottom  drag 
coefficient  similar  to  the  Manning  bottom  drag  used  by  Guan  et  al.  (2013),  who  used  a  Manning 
coefficient  of  0.013  s/m1/3.  Both  the  NCOM  and  Manning  bottom  drag  formulations  give  a  bottom 
drag  coefficient  of  about  0.0038  for  a  water  depth  of  0.08  nr  and  a  bottom  drag  coefficient  of  about 
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Fig.  3  —  SSH  from  model  simulation  of  dam-break  flow  over  a  triangular  bump  at  5-s  intervals. 
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Fig.  4  —  Comparison  of  observed  (blue  dots)  and  model-simulated  (black  line)  SSH  at  the  7  measurement  locations 
shown  in  Fig.  2  for  the  dam-break  flow  over  a  triangular  bump. 
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Fig.  5  —  Setup  and  dimensions  for  laboratory  experiment  for  a  dam-break  flow  into  a  channel  and  around  a  90° 
corner.  Top  and  side  views  are  shown.  The  gate  between  the  tank  and  the  channel  is  drawn  in  dark  blue.  The  six 
SSH  measurement  locations  are  shown  in  purple. 
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Fig.  6  —  Comparison  of  observed  (blue  dots)  and  model-simulated  (black  line)  SSH  at  the  6  measurement  locations 

shown  in  Fig.  5  for  the  dam-break  flow  around  a  90°  corner. 


0.0045  for  a  water  depth  of  0.05  m.  The  minimum  water  depth  Dmin  for  the  WAD  was  set  to 
0.01  m.  The  time  step  used  was  0.002  s.  The  simulation  was  run  for  60  s  and  output  from  the 
simulation  was  saved  at  0.1-s  intervals. 

Figure  6  shows  a  comparison  of  the  model  simulated  SSH  with  the  measured  values  at  the  six 
SSH  measurement  locations.  The  measured  values  correspond  to  the  laboratory  experiment  that 
was  run  with  a  “wet”  bed,  with  an  initial  water  depth  of  0.01  m  within  the  channel.  The  results  at 
gauges  3-6  in  Fig.  6  are  fairly  similar  to  those  obtained  by  Guan  et  al.  (2013)  (the  latter  did  not 
show  results  for  gauges  1  and  2). 

The  SSH  at  gauge  1  shows  the  simulated  drop  of  the  water  level  within  the  tank  to  be  in  good 
agreement  with  the  observed  drop.  The  SSH  at  gauges  2-4  within  the  first  section  of  the  channel 
shows  the  time  of  arrival  of  the  wetting  front  at  each  of  these  locations,  and  also  a  later  jump  in 
the  SSH  due  to  the  arrival  of  an  upstream  propagating  wave  caused  by  reflection  of  the  original 
wetting  front  at  the  end  of  the  first  section  of  the  channel.  The  time  of  arrival  of  the  reflected 
wave  at  gauges  4,  3,  and  2  in  the  first  section  of  the  channel  is  about  6,  11,  and  15  s,  respectively. 
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The  worst  agreement  between  the  simulated  and  measured  SSH  in  Fig.  6  is  at  gauge  2  up  until 
the  arrival  of  the  reflected  wave  at  about  15  s;  after  that  the  simulated  SSH  agrees  well  with  the 
measured  value.  A  similar  discrepancy  between  the  simulated  and  measured  SSH  at  gauge  2  was 
noted  by  Soares  Frazao  and  Alcrudo  (1999),  and  can  be  seen  in  Soares  Frazao  and  Zech  (1999)  in 
their  Fig.  6b.  The  SSH  at  gauges  5  and  6  in  the  second  section  of  the  channel  shows  the  arrival  of 
the  original  wetting  front  and  the  gradual  drop  of  the  SSH  as  the  water  drains  away. 

The  model  results  were  found  to  be  fairly  sensitive  to  the  bottom  drag,  which  was  noted  by 
Soares  Frazao  and  Alcrudo  (1999).  With  a  lower  bottom  drag,  the  flow  of  water  down  the  channel 
increases  in  speed,  the  water  level  within  the  tank  drops  more  quickly,  the  upstream  propagation 
of  the  wave  reflected  at  the  end  of  the  first  section  of  the  channel  is  slowed  because  of  the  increased 
speed  of  the  flow  in  the  channel,  and  the  water  drains  out  of  the  second  section  of  the  channel  more 
quickly.  The  latter  effect  is  especially  noticeable  in  lower  water  levels  at  gauge  6  after  the  arrival 
of  the  wetting  front.  All  of  these  changes  from  the  use  of  a  lower  bottom  drag  tend  to  reduce  the 
agreement  with  the  measured  SSH  relative  to  the  results  in  Fig.  6. 

3.5  San  Francisco  Bay 

Most  of  San  Francisco  Bay  (SFB)  is  relatively  shallow,  and  the  bay  has  a  variety  of  areas  that 
are  subject  to  WAD.  For  these  reasons  and  the  fact  that  we  had  done  some  earlier  work  in  SFB 
during  which  we  had  realized  that  a  WAD  capability  was  needed  to  conduct  proper  simulations  in 
this  area,  SFB  was  used  as  the  primary,  realistic,  test  case  for  testing  during  the  development  of 
NCOM’s  WAD  scheme.  Simulations  of  SFB  were  run  at  resolutions  of  500,  200,  and  100  m  and 
in  single-layer  barotropic  mode  and  in  multi-layer  baroclinic  mode.  The  coarser  resolution  runs 
allowed  for  quicker  turn-around,  while  the  higher  resolution  runs  provided  better  resolution  of  the 
finer-scale  features  of  the  bay  and  broader,  more  extensive,  WAD  areas  in  terms  of  the  number  of 
grid  points  involved. 

Note  that  there  are  a  number  of  low-lying  areas  around  SFB  that  are  not  regularly  flooded 
at  high  tide  because  they  are  protected  by  dikes.  Many  of  these  areas  were  once  part  of  the  bay, 
e.g.,  marshlands,  but  have  since  been  reclaimed  for  various  uses.  Since  the  dikes  are  often  fairly 
narrow,  they  usually  aren’t  well  resolved  in  bathymetric  data  bases  for  SFB.  Hence,  if  one  computes 
a  bathymetry  for  SFB  from  such  data  bases,  the  WAD  areas  may  appear  to  be  more  extensive  than 
they  actually  are.  One  can  examine  some  of  the  low-lying  areas  derived  from  the  bathymetric 
data  bases  using,  for  example,  Google  Maps,  and  see  there  is  human  activity  occurring  in  many  of 
these  areas,  e.g.,  farms,  buildings,  etc.,  that  would  not  be  there  if  the  areas  were  regularly  flooded. 
However,  for  the  purpose  of  testing  the  ability  of  an  ocean  model  to  perform  WAD  in  a  variety 
of  situations,  whether  an  area  does  actually  flood  at  high  tide  may  not  be  all  that  relevant,  i.e. , 
having  more  extensive  WAD  areas  tends  to  provide  a  better  test  of  the  robustness  of  the  WAD 
ability  of  the  model,  except,  of  course,  if  comparing  the  model  results  with  observations  in  or  near 
these  areas. 

Figure  7  shows  the  domain  and  bathymetry  used  for  most  of  the  SFB  simulations.  The  depths 
are  with  respect  to  mean  sea  level  (MSL).  The  two  plots  in  Fig.  7  show  the  depths  contoured 
between  zero  and  -100  m,  and  between  -8  and  +2m,  the  latter  to  better  illustrate  the  bathymetry 
in  the  shallow  areas  of  the  bay.  The  maximum  tidal  amplitudes  in  SFB  are  about  1  m.  Grid  cells 
with  an  elevation  greater  than  +2  m  are  defined  to  be  land  areas  that  are  outside  the  computational 
domain  and  these  areas  are  shown  in  black. 
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Fig.  7  —  Domain  and  bathymetry  used  for  SFB  simulations.  The  plot  on  the  left  shows  the  depths  contoured  between 
zero  and  -100  m  depth.  The  locations  of  IHO  stations  are  labeled.  The  plot  on  the  right  shows  the  depths  contoured 
between  -8  and  +2  m  to  better  illustrate  the  bathymetry  in  the  shallow  areas. 


SFB  is  composed  of  three  main  areas:  the  northwestern  part  is  San  Pablo  Bay,  the  northeastern 
part  is  Suisun  Bay,  and  the  southern  part  is  sometimes  referred  to  as  South  Bay.  The  San  Pablo 
and  Suisun  Bays  are  connected  by  the  Carquinez  Strait.  The  greatest  depths  within  SFB  are  about 
-100  nr  in  the  mouth  of  the  bay,  where  there  is  extensive  bottom  scouring  due  to  the  strong  tidal 
currents.  Most  of  the  rest  of  SFB  outside  the  main  channels  is  fairly  shallow,  i.e.,  less  than  5  nr  deep. 
As  noted  previously,  some  of  the  low-lying  areas  are  protected  from  flooding  by  dikes.  However, 
some  of  these  protected  areas  are  not  accounted  for  in  Fig.  7,  since  the  dikes  are  not  very  well 
resolved  in  the  bathymetric  data  base  that  was  used,  and  having  more  extensive  WAD  areas  was 
considered  to  be  useful  for  the  purpose  of  testing  the  WAD.  Another  issue  with  the  bathymetric 
data  base  used  for  Fig.  7  is  that  the  precision  of  the  land  elevations  is  only  to  the  nearest  meter. 
This  can  be  seen  in  Fig.  7b  by  the  lack  of  shading  of  the  color  contours  in  the  areas  that  have  an 
elevation  above  0  m. 

To  generate  fluctuating  water  levels  within  SFB,  tidal  forcing  was  applied  at  the  open  bound¬ 
aries,  which  lie  outside  SFB,  for  the  eight  main  tidal  constituents:  Ki,  Oi,  Pi,  Qi,  K2,  Mo,  N2, 
and  S2.  This  tidal  forcing  (i.e.,  the  tidal  SSH  and  current  transports)  was  obtained  from  one  of 
the  tidal  data  bases  (i.e.,  the  1/12°  West  Coast  data  base)  developed  at  Oregon  State  University 
(OSU)  by  Egbert  and  Erofeeva  (2003).  Tidal  potential  forcing  within  the  interior  of  the  domain 
was  found  to  have  little  effect  on  the  tides  in  the  simulations  because  of  the  relatively  small  size  of 
the  domain. 

During  the  testing  of  NCOM’s  WAD  scheme,  over  140  simulations  in  SFB  were  conducted 
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looking  at  different  bathymetries,  different  horizontal  grid  resolution  (500,  200,  and  100  m),  different 
vertical  grid  resolution  (single  and  multiple  layers),  different  types  of  vertical  grids  (hybrid  sigma/z- 
level  and  generalized  vertical  coordinate),  different  physics  (barotropic  and  baroclinic),  different 
input  parameters,  different  values  of  the  minimum  allowed  water  depth  Dmin  for  WAD  (10,  5,  2,  1, 
and  0.1  cm),  different  numbers  of  computer  processors,  and  different  numerical  precision  (single  and 
double  precision).  All  the  simulations  ran  smoothly  as  long  as  (i)  the  time  step  was  sufficiently  small 
to  prevent  instability  due  to  exceeding  the  advective  CFL  limits,  and  (ii)  the  errors  we  found,  which 
were  related  to  the  changes  being  made  for  the  WAD,  were  corrected  (many  of  these  simulations 
were  conducted  during  the  development  of  the  WAD  scheme). 

Figure  8  shows  a  comparison  of  the  model-simulated  tide  with  the  observed  tide  at  the  Inter¬ 
national  Hydrographic  Office  (IHO)  stations  within  SFB  shown  in  Fig.  7.  The  model-simulated 
tides  are  for  a  40-layer,  fully  baroclinic  simulation  on  the  200-m  SFB  grid.  The  comparison  is  fairly 
good  except  at  Oakland  Airport,  which  in  the  model  simulation  is  in  an  area  of  severely  restricted 
circulation  due  to  the  presence  of  some  islands  and  very  shallow  depths  (Fig.  7).  If  the  location  of 
the  model  point  used  to  represent  the  Oakland  Airport  IHO  station  is  moved  southwest  out  of  the 
area  of  restricted  flow,  the  model-simulated  tide  agrees  well  with  the  observed  tide.  It  is  likely  that 
the  bathymetry  near  this  IHO  station  is  not  adequately  resolved  and/or  is  too  shallow. 

Also,  the  model-simulated  tide  at  the  three  IHO  stations  in  Suisun  Bay  in  the  northeastern 
part  of  SFB  do  not  agree  with  the  observed  tide  as  well  as  at  most  of  the  other  stations.  This  may, 
in  part,  be  due  to  the  fact  that  the  model  domain  (a)  over-represents  the  WAD  areas,  and  (b)  has 
a  closed  boundary  on  the  eastern  end  of  Suisun  Bay  and  does  not  include  the  extension  of  the  delta 
areas  of  the  Sacramento  and  San  Joaquin  Rivers  further  to  the  east. 

The  total  area  of  the  SFB  domain  in  Fig.  7  that  is  within  the  computational  domain,  i.e. ,  not 
including  the  black-colored  areas  in  Fig.  7,  is  about  2820  km2.  For  the  tidal  simulations,  the  area 
that  always  remains  wet  is  about  2160  km2,  the  area  that  always  remains  dry  is  about  270  km2, 
and  the  area  of  WAD,  i.e.,  the  area  that  is  sometimes  wet  and  sometimes  dry,  is  about  390  km2. 
These  areas  represent  76%,  10%,  and  14%  of  the  computational  domain,  respectively. 

Figure  9  shows  sea-surface  temperature  (SST)  and  sea-surface  salinity  (SSS)  within  SFB  from 
a  500-m-resolution,  multi-layer,  baroclinic  simulation  with  variable  temperature  and  salinity.  The 
temperature  and  salinity  were  initialized  from  the  climatological  annual  mean  values.  The  forcing 
for  this  simulation  is  from  the  tides  and  river  inflows.  The  river  inflows  are  the  annual  mean  for 
the  Sacramento  and  San  Joaquin  Rivers,  which  are  about  682  and  213  m3/s,  respectively.  These 
are,  by  far,  the  two  largest  rivers  flowing  into  SFB  and  both  flow  into  the  eastern  end  of  Suisun 
Bay.  The  bathymetry  used  for  this  simulation  is  different  from  that  in  Fig.  7  in  that  the  areas  cut 
off  from  tidal  flooding  by  dikes  are  mostly  accounted  for;  hence,  the  horizontal  extent  of  the  WAD 
areas  is  significantly  reduced. 

The  SST  in  Fig.  9  is  at  10  d  into  the  simulation  and  shows  the  cooling  of  the  water  in  the  mouth 
of  SFB  due  to  vertical  mixing  by  the  strong  tidal  currents,  which  at  this  time  are  ebbing.  The  SSS 
in  Fig.  9  is  at  100  d,  which  is  when  the  salinity  distribution  within  the  bay  in  this  simulation  has 
almost  reached  equilibrium.  The  salinity  in  Suisun  Bay  in  Fig.  9  is  fairly  low  due  to  the  large  river 
inflows. 

A  summary  of  some  of  the  results  from  the  simulations  in  SFB  with  respect  to  the  WAD  are 
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Fig.  8  —  Comparison  of  observed  (solid  line)  and  model-simulated  (dashed  line)  tide  at  the  location  of  several  IHO 
tidal  stations  within  SFB.  The  time  period  of  the  plots  is  January  6-10,  2004. 
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Fig.  9  —  SST  in  °C  at  10  d  (left)  and  SSS  in  psu  at  100  d  (right)  for  a  multi-layer,  baroclinic  simulation  of  SFB  with 

river  inflows. 


as  follows.  Simulations  tend  to  run  more  smoothly  with  larger  time  steps  when  the  bathymetry 
is  smoother  due  to  lower  values  of  the  maximum  advective  CFL  values.  Baroclinic  simulations 
with  multiple  layers  tend  to  require  smaller  time  steps  as  the  number  of  layers  is  increased  due  to 
the  decreasing  thickness  of  the  layers  and  the  increased  possibility  of  exceeding  the  CFL  limit  for 
vertical  advection.  Double  precision  (i.e. ,  more  than  32-bit  precision  for  real  numbers)  is  needed  for 
accurate  conservation  of  scalar  fields  on  our  workstations,  especially  for  long  simulations  without 
data  assimilation,  and  for  high-resolution  simulations  with  a  small  time  step  (since  the  cumulative 
effects  of  roundoff  error  tend  to  increase  with  the  number  of  time  steps).  Both  the  advection  and 
vertical  mixing  terms  tend  to  be  a  source  of  roundoff  error,  i.e.,  fixing  just  one  or  the  other  does 
not  solve  the  roundoff  error  problem.  If  a  scalar  field  (e.g.,  salinity)  is  initialized  to  a  constant 
value,  and  the  values  transported  in  or  out  through  the  boundaries  and  through  source  or  sink 
terms  in  the  interior  are  the  same,  constant  value,  then  the  constant  value  of  the  scalar  field  will  be 
maintained  to  good  accuracy  if  double  precision  is  used  (this  is  a  useful  test  for  tracer  conservation) . 
The  use  of  the  Flux-Corrected- Transport  (FCT)  scheme  for  the  advection  of  scalar  fields  is  effective 
for  preventing  advective  overshoots.  The  sigma/z- level  (SIGZ)  and  generalized  vertical  coordinate 
(GVC)  versions  of  the  NCOM  code  give  almost  identical  results  for  the  same  model  setup  (note 
that  identical  results  are  not  possible  due  to  unavoidable  differences  in  the  sequence  of  some  of  the 
calculations) . 

3.6  Chesapeake  Bay 

Chesapeake  Bay  (ChB)  is  referred  to  as  a  “partially  mixed”  estuary.  Relatively  high  salinity 
water  from  the  ocean  flows  into  and  northward,  up  the  bay,  near  the  bottom,  and  fresher  water 
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from  a  number  of  rivers  flows  into  the  bay  and  southward,  down  the  bay,  near  the  surface  and 
out  of  the  mouth  of  the  bay.  This  type  of  circulation  within  a  bay  is  referred  to  as  an  “estuarine” 
circulation. 

Due  to  mixing  between  the  inflowing,  higher  salinity  water  near  the  bottom  and  the  fresher 
water  near  the  surface,  there  is  a  gradual  decrease  in  the  salinity  of  the  near-bottom,  saltier  water 
as  it  flows  up  the  bay  and  a  gradual  increase  in  the  salinity  of  the  fresher  water  as  it  flows  southward 
near  the  surface  towards  the  mouth  of  the  bay.  In  spite  of  this  vertical  mixing,  at  most  locations 
within  ChB,  there  is  a  significant  vertical  salinity  stratification  (Stroup  and  Lynn  1963);  hence,  the 
estuary  is  referred  to  as  being  “partially  mixed” . 

The  annual  mean  freshwater  inflow  from  the  rivers  into  ChB  is  about  2000  m3/s  and  the  mean 
of  the  inflowing  seawater  into  the  bay  is  on  the  order  of  7000  m3/s  (Valle-Levinson  et  al.  1998), 
so  that  the  mean  of  the  outflow  from  the  bay  is  on  the  order  of  9000  m3/s.  Simulations  of  the 
estuarine  circulation  in  ChB  with  NCOM  by  Martin  (1999)  developed  estuarine  flow  volumes  and 
a  salinity  structure  within  ChB  similar  to  the  observed  values.  Note  that  the  salinity  structure  in 
this  simulation  took  over  a  year  to  reach  a  fairly  steady  state  (Martin  1999). 

Recently,  numerical  model  simulations  were  conducted  by  the  Naval  Research  Laboratory 
(NRL)  Code  7320  in  the  Chesapeake  Bay  area  in  support  of  the  Trident  Warrior  2013  (TW13) 
Naval  Exercise  (Allard  et  al.  2014).  Since  the  tidal  amplitude  in  this  area  tends  to  be  less  than 
1  nr,  it  was  decided  to  set  the  shallowest  bottom  depth  H  at  sea  points  to  -2  m  for  the  simulations 
to  avoid  drying  out  grid  cells  during  the  ocean  model  runs.  This  limit  worked  satisfactorily  during 
the  exercise  itself,  which  took  place  during  the  summer  when  the  winds  tend  to  be  light.  However, 
during  spin  up  of  the  ocean  model  domain,  which  was  begun  in  January  2013,  strong  winds  during 
the  winter  caused  drying  out  of  some  areas  within  ChB  and  along  the  coast  due  to  set  down  of  the 
SSH  in  excess  of  2  m. 

Since  a  beta-test  version  of  NCOM  with  WAD  was  available  at  this  time,  it  was  decided  to  use 
the  new  WAD  capability  of  NCOM  for  the  spin  up  to  avoid  having  to  deepen  the  shallowest  allowed 
bottom  depth  beyond  the  -2  m  that  had  originally  been  used  to  set  up  the  ocean  model  domains  for 
TW13.  Hence,  NCOM  with  WAD  was  used  to  spin  up  the  ocean  circulation,  and  this  spinup  was 
then  used  to  initialize  some  of  the  TW13  ocean  forecasts.  The  spinup  was  run  from  the  beginning 
of  January  to  the  end  of  June  2013.  Note  that,  as  mentioned  above,  this  is  too  short  a  time  to  allow 
full  development  of  the  salinity  structure  within  ChB;  hence,  the  salinity  structure  within  and  just 
outside  ChB  for  this  spinup  was  initialized  from  a  2-year  spinup  from  earlier  work  with  NCOM  in 
ChB.  This  was  blended  with  salinity  values  in  the  deeper  water  from  Global  HYCOM,  which  also 
provided  the  initial  conditions  for  the  SSH,  ocean  current,  and  temperature  fields. 

Forcing  for  the  spinup  consisted  of  the  following:  monthly  climatological  river  inflows  for  the 
seven  largest  rivers  flowing  into  ChB  (the  James,  York,  Rappahannock,  Potomac,  Patauxent,  Pat- 
apsco,  and  Susquehanna)  with  a  combined  annual  mean  inflow  of  2220  m3/s  and  the  Delaware  River 
in  Delaware  Bay  with  an  annual  mean  inflow  of  570  m3/s,  tidal  forcing  at  the  open  boundaries  for 
the  eight  main  tidal  constituents  (Ki,  Oi,  Pi,  Qi,  K2,  M2,  N2,  S2)  from  the  OSU  1/12°  Atlantic 
Ocean  tidal  database  and  tidal  potential  forcing  over  the  interior  for  the  same  eight  constituents, 
3-hourly  values  of  atmospheric  forcing  from  the  operational  COAMPS  West  Atlantic  model  (con¬ 
sisting  of  surface  atmospheric  pressure,  wind  stress,  and  solar  and  net  longwave  radiation,  and 
latent  and  sensible  surface  heat  fluxes  computed  using  the  COAMPS  winds  and  air  temperature 
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and  humidity  and  the  NCOM  SST,  which  provides  some  feedback  to  the  surface  heat  fluxes  from 
the  NCOM  SST),  and  daily  boundary  conditions  for  the  SSH,  current  velocity,  temperature,  and 
salinity  fields  from  Global  HYCOM. 

Figure  10  shows  the  domain  used  for  the  spin  up.  The  grid  is  115  x  226  points  and  the  horizontal 
resolution  is  2  km.  The  relatively  coarse  grid  is  the  reason  these  plots  look  a  bit  blocky.  The  vertical 
grid  uses  a  total  of  40  layers,  with  19  sigma  layers  between  the  surface  and  the  bottom  down  to  a 
maximum  depth  of  -105  m,  and  with  21  fixed-depth  levels  (z-levels)  between  -105  m  and  -2800  m. 
The  vertical  grid  was  smoothly  stretched  downward  from  the  surface,  with  an  upper-layer  thickness 
of  0.95%  of  the  total  depth  of  the  sigma  part  of  the  grid;  hence,  the  thickness  of  the  upper  layer  is 
about  1  m  in  water  deeper  than  -105  m  and  is  proportionally  less  in  water  shallower  than  -105m. 
Note  that  the  depths  in  Fig.  10a  are  only  contoured  between  zero  and  -100  m  to  provide  more  detail 
in  the  bays  and  on  the  shelf.  The  depths  in  this  domain  increase  substantially  past  the  edge  of  the 
continental  shelf.  The  maximum  depths  within  this  domain  are  about  -2800  m  near  the  southeast 
corner. 

Figure  10b  shows  water  depths  within  the  computational  domain  contoured  within  the  range 
of  -4  to  -2  m  to  illustrate  areas  of  potential  WAD.  Without  WAD,  set  down  of  the  SSH  and  drying 
out  of  grid  cells  in  some  of  these  shallow  areas  caused  the  model  to  crash.  With  WAD,  the  model 
spinup  ran  smoothly  through  the  entire  January-June  period  without  a  problem  each  of  the  several 
times  it  was  run  (several  different  spinups  were  conducted  with  different  bathymetries) . 

Figure  11  shows  the  SST  and  SSS  at  the  end  of  the  spinup  of  the  ChB  domain  on  30  June  2013. 
The  SST  shows  some  cooling  along  the  coast  due  to  the  moderately  strong  south  winds  that  were 
occurring  at  this  time.  The  SSS  shows  the  large  variations  of  salinity  in  ChB  and  Delaware  Bay. 
The  SSS  at  the  head  of  these  bays  is  fairly  low  and  increases  towards  the  mouth  of  the  bays  due  to 
mixing  with  saltier  water  traveling  up  the  bays  near  the  bottom.  The  south  winds  and  upwelling 
along  the  coast  have  reduced  the  SSS  signature  of  the  outflow  plume  from  ChB  outside  the  mouth 
of  the  bay. 

We  refer  to  the  setup  used  here  for  ChB  as  an  example  of  “occasional”  WAD,  i.e.,  the  focus 
was  not  on  simulating  the  WAD  of  extensive  areas,  but  of  maintaining  robust  performance  of 
the  ocean  model  during  the  occasional  times  when  relatively  shallow  grid  cells  dry  out.  The 
purpose  of  including  this  simulation  of  ChB  in  the  VTR  was  to  illustrate  the  use  of  such  occasional 
WAD  to  maintain  numerical  stability  in  a  relatively  long  (6-month)  NCOM  simulation.  Since  the 
WAD  procedure  was  only  invoked  occasionally,  there  was  fairly  negligible  additional  computer  time 
required  for  the  WAD  calculations  that  were  needed. 

3.7  Cook  Inlet,  Alaska 

Cook  Inlet  (Cl)  is  located  on  Alaska’s  south-central  coast.  The  upper  part  of  Cl  has  a  tidal 
range  of  over  10  m,  and  these  are  the  second  highest  tides  in  North  America  after  the  Bay  of  Fundy 
in  Nova  Scotia.  In  part  due  to  the  large  tides,  Cl  has  extensive  WAD  areas,  which  are  composed 
mostly  of  tidal  mud  flats  that  are  exposed  during  low  tide.  Cl  has  been  used  for  testing  both 
tidal  prediction  and  WAD  in  numerical  models  by  other  investigators,  e.g.,  Matthews  and  Mungall 
(1972),  Oey  et  al.  (2007),  and  Kowalik  and  Proshutinsky  (2010). 

Figure  12  shows  the  domain  and  the  bathymetry  used  for  the  numerical  simulations  of  the  tides 
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Fig.  10  -  Domain  and  bathymetry  used  for  Chesapeake  Bay  simulations.  The  plot  on  the  left  shows  the  depths 
contoured  between  zero  and  -100  m  depth.  The  plot  on  the  right  shows  the  depths  contoured  just  between  -2  and  -4 
m  to  indicate  areas  where  WAD  can  potentially  occur. 
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Fig.  11  —  SST  (left)  in  °C  and  SSS  (right)  in  psu  at  the  end  of  the  spinup  of  the  ChB  domain  on  30  June  2013 
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in  CL  The  town  of  Anchorage,  Alaska  is  located  in  the  northeastern  part  of  Cl  where  the  waters  of 
Cl  separate  into  Knik  Arm  to  the  northeast  and  Turnagain  Arm  (TA)  to  the  east-southeast.  Note 
that  the  eastern  two  thirds  of  TA  is  relatively  shallow  and  at  low  tide  consists  mostly  of  exposed 
tidal  mudflats.  As  the  rising  tide  enters  the  eastern  two  thirds  of  TA,  a  tidal  bore  is  usually 
generated  by  the  rapidly  rising  tide,  which  propagates  eastward  up  TA.  The  island  just  west  of 
Anchorage  is  called  Fire  Island,  and  at  low  tide  it  is  possible  to  walk  from  Anchorage  to  Fire  Island 
on  the  exposed  mud  flats,  though  the  rapid  advance  of  the  flooding  tide  and  the  quicksand-like 
properties  of  occasional  soft  areas  on  the  mudflats  can  be  treacherous  both  here  and  in  most  of  the 
other  WAD  areas  in  CL 

The  bathymetry  for  Cl  in  Fig.  12  was,  for  the  most  part,  obtained  from  the  National  Oceanic 
and  Atmospheric  Administration  (NO  A  A)  Digital  Environmental  Model  (DEM)  Number  937, 
which  has  a  resolution  of  24  s  (about  1  km)  and  is  one  of  a  number  of  such  data  bases  being 
developed  by  NOAA  for  tsunami  and  flood  prediction  and  is  available  online.  However,  this  and 
other  bathymetry  data  bases  we  investigated  for  Cl  had  poor  accuracy  in  many  areas,  notably  in 
the  WAD  areas  and  in  the  NE  part  of  Cl  including  Knik  Arm  and  TA.  This  is  and  has  been  a 
significant  problem  for  those  trying  to  model  the  tides  and  WAD  in  CI. 

However,  we  were  able  to  obtain  some  more  accurate,  high-  (50-100  m)  resolution,  bathymetry 
data  for  the  NE  part  of  CI  from  Tal  Ezer  (personal  communication)  and  these  data  are  included  in 
Fig.  12.  Tal  Ezer  and  Hua  Liu  have  been  working  on  developing  a  bathymetry  for  this  area  using 
a  combination  of  satellite  photographs  and  tidal  data  to  determine  the  location  of  the  land-sea 
boundary  at  various  phases  of  the  tide  and  thereby  deduce  the  bathymetry  (Ezer  and  Liu  2009, 
2010).  The  bathymetry  data  developed  by  Ezer  and  Liu  (EL)  includes  depths  in  the  WAD  areas 
estimated  from  the  satellite  and  tidal  data  and  in  the  deeper  areas  from  nautical  charts. 

The  horizontal  grid  resolution  used  for  the  tidal  simulations  in  CI  was  about  1  km  and  the 
horizontal  dimensions  of  the  longitude-latitude  grid  were  312  x  356  points.  Only  a  single  layer 
was  used  in  the  vertical.  A  quadratic  form  of  the  bottom  drag  was  used  (see  Section  2.5)  with  a 
minimum  value  for  the  bottom  drag  coefficient  of  0.0025  and  a  bottom  roughness  of  0.01  m.  The 
third-order  upwind  scheme  was  used  for  momentum  advection.  Tidal  forcing  was  from  the  OSU 
1/12°  Pacific  Ocean  data  base.  The  eight  main  tidal  constituents  (Ki,  Oi,  Pi,  Qi,  K2,  M2,  N2, 
and  S2)  were  used  for  the  tidal  forcing.  The  time  step  used  was  30  s. 

Figure  13  shows  a  comparison  of  the  predicted  tide  with  the  tide  computed  at  eight  IHO  tidal 
stations  within  and  just  outside  CI  (the  locations  of  these  IHO  stations  are  noted  in  Fig.  12).  The 
observed  tide  was  computed  using  the  IHO  data  for  the  same  eight  tidal  constituents  that  were 
used  for  the  model  simulation.  The  time  period  in  Fig.  13  is  the  beginning  of  day  6  through  day 
10  of  the  model  simulation  (i.e.,  January  6  through  10,  2001).  (All  the  times  referred  to  here  are 
GMT.)  Note  that  this  time  period  includes  a  sampling  of  the  highest  tides  regularly  observed  in 
CI.  The  close  agreement  at  the  Seward  tidal  station  in  Fig.  13d,  which  is  outside  CI  on  the  south 
facing  coast  east  of  the  mouth  of  CI  (Fig.  12),  suggests  that  the  accuracy  provided  by  the  OSU 
Pacific  tidal  database  that  was  used  to  provide  the  tidal  boundary  conditions  is  very  good. 

The  overall  agreement  of  the  observed  and  simulated  tides  in  Fig.  13  is  quite  good.  There  is 
a  slight  lag  in  the  predicted  tidal  phase  and  a  slight  overprediction  of  the  tidal  amplitude  at  the 
Fire  Island  and  Anchorage  tidal  stations  in  the  NE  part  of  CI.  Both  of  these  errors  could  probably 
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Fig.  12  Bathymetry  for  Cook  Inlet.  Depths  are  relative  to  mean  sea  level.  Areas  with  land  elevations  greater  than 
8  m  are  outside  the  computational  domain  and  are  shown  in  black.  The  locations  of  eight  IHO  stations  are  labeled. 
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Fig.  13  —  Comparison  of  observed  (solid  line)  and  model-simulated  (dashed  line)  tide  at  several  IHO  stations  within 
Cook  Inlet.  The  time  period  of  the  simulations  is  January  6  through  10,  2001. 
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be  reduced  by  increasing  the  damping  of  the  tide  as  it  propagates  up  Cl  by  slightly  increasing  the 
bottom  drag. 

There  are  many  WAD  areas  within  Cook  Inlet  and  along  the  coast  just  outside  of  Cook  Inlet. 
As  noted  earlier,  the  bathymetry  in  these  areas  is,  with  respect  to  the  accuracy  needed  for  validating 
a  WAD  model,  not  very  well  known.  However,  within  the  narrow,  upper  two-thirds  of  Turnagain 
Arm  (TA)  in  northeast  Cl,  the  bottom  is  almost  fully  exposed  at  low  tide  and  there  is  a  tidal  bore 
that  propagates  up  TA  on  the  flooding  tide.  The  tidal  bore  is  a  popular  attraction  in  this  part 
of  Alaska,  and  can  easily  be  observed  from  the  Seward  Highway  that  runs  along  the  north  side  of 
TA.  The  approximate  times  of  occurrence  of  the  bore  at  several  locations  along  TA  are  provided 
to  assist  those  wishing  to  view  the  bore.  Because  of  the  availability  of  this  information  on  the  tidal 
flooding  that  occurs  in  TA,  it  was  decided  to  focus  the  validation  of  WAD  in  Cl  on  the  WAD  that 
occurs  in  TA. 

Figure  14  shows  the  bathymetry  in  the  northeast  part  of  Cl,  which  includes  both  Knik  Arm 
and  TA,  provided  by  Tal  Ezer.  The  bathymetry  is  contoured  between  -5  and  +5  m  to  highlight 
the  WAD  areas.  Figure  14  shows  the  locations  of  the  Fire  Island  and  Anchorage  IHO  stations  and 
the  names  of  several  landmarks  along  the  northern  side  of  TA. 

The  time  of  passage  of  the  tidal  bore  as  it  propagates  up  TA  is  available  at  the  locations  along 
TA  labeled  in  Fig.  14  relative  to  the  time  of  low  tide  at  Anchorage.  These  observed  times  are  listed 
in  Table  1.  For  the  model  simulations,  the  times  of  arrival  of  the  tidal  bore  listed  in  Table  1  are 
taken  to  be  the  time  at  which  the  SSH  exceeds  the  minimum  SSH  during  the  previous  low  tide  at 
that  location  by  0.5  m.  Note  that  the  values  in  Table  1  for  the  model  simulations  are  an  average 
over  nine  low  tides  between  January  6  and  10,  2001. 

Table  1  -  Comparison  of  observed  and  model-predicted 
tidal  bore  arrival  times  in  Turnagain  Arm  with  respect 
to  time  of  low  tide  at  Anchorage.  Values  are  in  minutes 
and  are  an  average  over  nine  low  tides  between  January 
5  and  10,  2001. 
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Calculation  of  the  arrival  times  for  the  tidal  bore  for  our  initial  tidal  simulation  for  Cl  found 
the  simulated  bore  to  lag  far  behind  the  observed  bore  (see  Expt.  1  in  Table  1).  For  this  simulation, 
the  tidal  bore  did  not  propagate  very  far  beyond  Indian  Pt.  and  so  no  arrival  times  are  given  for 
locations  beyond  Indian  Pt.  for  Expt.  1  in  Table  1. 


In  order  to  try  to  increase  both  the  speed  and  extent  of  the  tidal  flooding  in  TA,  the  bottom 
drag  was  significantly  reduced,  i.e,  the  minimum  value  of  the  bottom  drag  coefficient  Cbmin  was 


31 


Latitude  (N) 


Martin  et  al. 


61.53 


5 

4 

3 

2 

1 

0 

-1 

-2 

-3 

-4 

-5 

i  rdwood 


Longitude  (E) 


60.81  -f 
209.6 


Fig.  14  —  Bathymetry  for  the  northeast  part  of  Cook  Inlet  provided  by  Tal  Ezer.  Depths  are  with  respect  to  MSL. 
The  locations  of  the  Fire  Island  and  Anchorage  IHO  stations  are  indicated  and  the  locations  of  several  landmarks 
along  the  northern  side  of  Turnagain  Arm  are  also  shown. 
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reduced  from  0.0025  to  0.0001  and,  just  within  the  shallow  areas  of  TA,  the  bottom  roughness  was 
reduced  from  0.01  to  the  very  small  value  of  10-20  nr  so  as  to  maintain  Cb  =  Cbmin  in  the  WAD 
areas  (see  Section  2.5).  This  reduction  of  the  bottom  drag  significantly  increased  the  extent  of  the 
tidal  flooding  in  TA,  but  still  resulted  in  very  large  lags  in  the  model-simulated  bore  times  relative 
to  the  observed  times  (see  Expt.  2  in  Table  1). 

One  possible  factor  contributing  to  the  slow  arrival  of  the  tidal  bore  in  TA  relative  to  the 
time  of  low  tide  at  Anchorage  is  the  presence  of  some  shallow  features  in  the  western  third  of  TA, 
including  a  shallow  ridge  connecting  the  mud  flats  on  the  north  and  south  sides  of  TA,  which  can 
be  seen  in  Fig.  14.  This  shallow  ridge  blocks  the  tidal  drainage  out  of  the  part  of  TA  east  of  the 
ridge,  and  also  delays  the  arrival  of  the  tidal  flood  east  of  the  ridge.  Hence,  the  bathymetry  in  the 
wider,  western  third  of  TA  was  modified  by  imposing  an  approximately  4-km  wide  channel  though 
the  middle  of  the  western  third  of  TA,  with  a  maximum  depth  of  -12  m  near  the  mouth  of  TA 
and  a  maximum  depth  of  -4  m  near  the  beginning  of  the  narrow,  eastern  two  thirds  of  TA.  These 
changes  in  the  bathymetry  allowed  for  improved  drainage  of  the  tide  and  a  quicker  arrival  of  the 
flood  just  east  of  the  (now  removed)  ridge.  However,  the  arrival  time  of  the  bore  at  locations  along 
the  narrow,  upper  two  thirds  of  TA  was  not  improved  (see  Expt.  3  in  Table  1).  This  delay  now 
appeared  to  be  due  mainly  to  the  large  increase  in  the  bottom  elevation  at  the  beginning  of  the 
narrow  part  of  TA  near  Beluga  Pt.,  which  rises  from  -4  m  to  0  m  over  a  distance  of  a  few  km. 

At  this  point,  it  was  decided  to  try  more  extensive  modification  of  the  bathymetry  in  the 
narrow,  eastern  two-thirds  of  TA.  The  minimum  depth  in  TA  near  Beluga  Pt.  in  Fig.  14  is  about 
0  nr  and  the  minimum  depth  rises  to  about  +2  m  near  Indian  Pt.  The  minimum  depth  of  0  m 
near  Beluga  Pt.  would  seem  to  make  it  difficult  for  the  tidal  flood  to  arrive  at  this  location  much 
before  the  tide  in  the  western  part  of  TA  rises  to  this  level,  which  would  be  about  3  hrs  after  the 
occurrence  of  low  tide.  This  analysis  is  consistent  with  the  arrival  times  of  the  bore  at  Beluga  Pt. 
for  Expts.  2  and  3  in  Table  1.  Another  problem  is  that  the  depth  in  the  narrow  part  of  TA  in 
Fig.  14  reaches  a  maximum  of  about  +3  m  at  Bird  Pt.  and  then  decreases  towards  the  end  of  TA 
near  Portage.  This  raises  the  question  of  how  the  area  above  Bird  Pt.  is  to  drain  if  the  bottom 
slopes  the  wrong  way. 

The  simplest  idealized  bathymetry  to  try  was  considered  to  be  a  linear  variation  of  the  bottom 
slope  along  the  entire,  narrow,  eastern  two-thirds  of  TA.  Hence,  a  linear  bottom  slope  was  tried 
with  a  minimum  depth  of  -4  m  at  the  beginning  of  the  narrow  part  of  TA  just  west  of  Beluga 
Pt.  and  a  minimum  depth  of  +3  m  at  the  end  of  TA  near  Portage.  This  gives  minimum  depths 
at  Beluga,  Indian,  and  Bird  Pts.  and  at  Girdwood  of  -3.4,  -2.9,  -0.3,  and  +1.3  nr,  respectively. 
Experiment  4  in  Table  1  shows  the  predicted  bore  times  in  TA  for  this  simulation.  The  results  with 
respect  to  the  earlier  experiments  are  noticeably  improved,  but  the  simulated  bore  still  significantly 
lags  the  observed  bore,  i.e.,  the  lags  at  Beluga,  Indian,  and  Bird  Pts.  and  at  Gridwood  are  38,  46, 
78,  and  43  minutes,  respectively.  Note  that  the  steady  increase  in  the  minimum  elevation  along 
the  upper  two-thirds  of  TA  in  this  simulation  does  allow  most  of  the  water  in  this  part  of  TA  to 
drain  out  at  low  tide. 

At  this  point  in  our  investigation  of  the  tidal  bore  in  TA,  Tal  Ezer  notified  us  that  a  new 
bathymetry  for  Cl  was  available  from  NOAA  (Zimmermann  and  Prescott  2014)  and  suggested  we 
look  at  it.  We  contacted  Mark  Zimmermann  of  NOAA’s  Alaska  Fisheries  Science  Center  in  Seattle, 
who  helped  us  obtain  this  bathymetry.  The  Zimmermann  and  Prescott  (ZP)  bathymetry  is  at  50-rn 
resolution  and  covers  almost  all  of  the  interior  of  Cl,  except  for  a  small  area  near  the  mouth  of 
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CL  The  ZP  bathymetry  was  interpolated  to  our  model  grid,  and  the  NOAA  DEM  937  bathymetry 
was  used  for  the  areas  in  our  domain  not  covered  by  the  ZP  bathymetry.  The  ZP  bathymetry  as 
received  was  referenced  to  mean  lower  low  water  (MLLW).  To  convert  the  vertical  datum  to  MSL, 
we  used  our  tidal  solution  for  Cl  along  with  the  difference  between  MSL  and  MLLW  at  several  tidal 
stations  within  Cl  provided  online  by  NOAA  (Table  2).  Note  that  the  correction  is  about  5  m  in 
the  NE  part  of  CI. 

Table  2  —  Difference  between  MSL  and  MLLW  (nr)  at 
several  IHO  tidal  stations  in  CI. 


Homer 

Seldovia 

Nikiski 

Anchorage 

1.69 

2.90 

3.44 

5.02 

Figure  15  shows  the  ZP  bathymetry  in  the  NE  part  of  CI  for  the  same  region  shown  in  Fig.  14 
for  the  bathymetry  from  Tal  Ezer.  Comparison  of  Fig.  15  with  Fig.  14  shows  significant  differences. 
For  the  ZP  bathymetry,  the  flow  through  the  western  third  of  TA  is  less  obstructed  by  shallow  areas, 
and  the  minimum  bottom  elevation  in  the  eastern  two-thirds  of  TA  up  to  Girdwood  is  significantly 
lower.  Both  of  these  differences  should  help  to  advance  the  arrival  time  of  the  tidal  bore  in  TA. 

With  this  revised  bathymetry  for  CI,  the  predicted  tidal  amplitudes  at  Fire  Island  and  An¬ 
chorage  were  about  0.8  nr  lower  than  observed  for  the  highest  tides.  Hence,  the  bottom  roughness 
in  the  main  part  of  Cook  Inlet  was  reduced  from  0.01  to  0.003  m,  which  resulted  in  increased  tidal 
amplitudes  in  the  NE  part  of  CI  that  agreed  better  with  the  observed  tides.  This  change  in  the 
bottom  roughness  in  the  main  part  of  CI  was  used  for  all  the  simulations  conducted  with  the  ZP 
bathymetry. 

Experiment  5  in  Table  1  shows  the  predicted  bore  times  for  a  tidal  simulation  with  the  ZP 
bathymetry.  The  time  of  the  bore  with  respect  to  low  tide  at  Anchorage  is  significantly  reduced 
relative  to  the  previous  experiments,  except  at  Portage  where  the  elevation  (4  nr)  is  now  higher 
than  in  the  previous  simulations.  However,  the  presence  of  some  isolated  deeper  pockets  in  the 
eastern  two  thirds  of  TA  (Fig.  15)  prevent  these  areas  from  draining  completely  during  the  tidal 
ebb.  The  water  depth  at  low  tide  in  some  of  these  pockets  exceeds  5  m. 

Figure  16  shows  plots  of  SSH  versus  time  for  Expt.  5  for  Anchorage  and  for  the  locations  along 
TA.  The  locations  along  TA  never  completely  dry  out  for  any  length  of  time,  except  at  Portage. 
The  SSH  at  Beluga,  Indian,  and  Bird  Pts.  shows  a  bit  of  noise  at  low  tide,  which  is  due  to  some 
residual  drainage  occurring  in  these  areas  between  the  remaining  pockets  of  water.  Because  of  the 
incomplete  drainage  in  Expt.  5,  some  simulations  were  conducted  with  the  ZP  bathymetry  with 
modified  depths  in  the  eastern  two-thirds  of  TA  to  try  to  achieve  more  complete  drainage. 

Experiment  6  in  Table  1  shows  the  predicted  bore  times  for  a  tidal  simulation  with  the  ZP 
bathymetry  with  the  minimum  depth  at  Beluga  Pt.  set  to  -5.5  m  and  minimum  depths  at  Indian 
and  Bird  Pts.  and  at  Girdwood  and  Portage  of  -4.55,  -1.5,  +0.5,  and  +5.1  m,  respectively.  These 
minimum  depths  yield  a  bottom  slope  along  the  thalweg  of  TA  between  Beluga  Pt.  and  Girdwood 
of  about  2  nr  per  10  km.  The  bore  times  are  generally  slightly  improved  over  Expt.  5. 

Figure  17  shows  plots  of  SSH  versus  time  at  Anchorage  and  along  TA  for  Expt.  6.  An  extensive 
dry  period  occurs  during  low  tide  at  Bird  Pt.,  Girdwood,  and  Portage,  as  indicated  by  the  flattening 
of  the  SSH  during  low  tide  at  these  locations.  At  Beluga  and  Indian  Pts.,  a  brief  dry  period  occurs 
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Fig.  15  -  Depths  in  the  northeast  part  of  Cl  from  the  Cl  bathymetry  of  Zimmermann  and  Prescott  (2014).  Depth 
are  relative  to  MSL.  The  locations  of  the  Fire  Island  and  Anchorage  IHO  stations  are  indicated  and  the  locations  of 
several  landmarks  along  the  north  side  of  Turnagain  Arm  are  also  shown. 
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Fig.  16  —  Plots  of  SSH  vs  time  for  Expt.  5  at  the  Anchorage  IHO  station  and  at  several  locations  along  TA  from 

January  6  through  January  10,  2001. 
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Fig.  17  —  Plots  of  SSH  vs  time  for  Expt.  6  at  the  Anchorage  IHO  station  and  at  several  locations  along  TA  from 

January  6  through  January  10,  2001. 

on  the  lowest  tides.  The  lack  of  noise  in  the  SSH  at  low  tide  compared  with  Fig.  16  indicates  that 
there  is  not  much  residual  drainage  occurring  at  low  tide.  The  increase  in  the  maximum  SSH  along 
TA  in  Fig.  17  indicates  the  effect  of  momentum  conservation  (inertia)  on  the  tidal  flood  in  TA.  This 
effect  is  partly  dependent  on  the  narrowing  of  TA  from  west  to  east.  The  increase  in  the  SSH  along 
TA  also  depends  on  the  bottom  drag  and  would  be  reduced  if  the  bottom  drag  were  increased. 

Experiment  7  in  Table  1  shows  the  predicted  bore  times  for  a  tidal  simulation  with  the  ZP 
bathymetry  with  the  minimum  depths  between  Beluga  Pt.  and  Girdwood  increased  by  0.5  m  over 
the  corresponding  depths  used  for  Expt.  6,  and  the  minimum  bottom  depth  at  Portage  increased 
by  0.2  m  to  5.3  m.  This  change  was  made  to  see  if  the  drainage  at  Beluga  Pt.  could  be  improved 
over  Expt.  6.  These  minimum  bottom  depths  increase  the  elevation  of  the  bottom  between  Beluga 
Pt.  and  Girdwood  by  0.5  m,  but  maintain  the  same  bottom  slope  in  this  area.  (It  was  found  that 
lower  bottom  slopes  in  this  area  tend  to  reduce  the  degree  of  drainage  at  Beluga  Pt.)  Table  1  shows 
the  bore  times  are  slightly  increased  by  this  change,  which  was  expected.  Plots  of  the  SSH  versus 
time  along  TA  (not  shown)  look  similar  to  the  results  for  Expt.  6  in  Fig.  17,  though  the  period  of 
drying  at  Beluga  Pt.  is  slightly  increased  over  Expt.  6. 
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The  total  area  of  the  Cl  domain  in  Fig.  12  that  is  within  the  computational  domain,  i.e. , 
not  including  the  black-colored  areas,  is  about  37200  km2.  For  the  Expt.  7  tidal  simulation,  the 
area  that  always  remains  wet  is  about  32350  km2,  the  area  that  always  remains  dry  is  about  1970 
km2,  and  the  area  of  WAD,  is  about  2880  km2.  These  areas  represent  87%,  5%,  and  8%  of  the 
computational  domain,  respectively. 

Looking  at  just  the  area  within  Cl  itself,  the  total  area  that  is  within  the  computational  domain 
is  about  21020  km2.  For  Expt.  7,  the  area  that  always  remains  wet  is  about  17180  km2,  the  area  that 
always  remains  dry  is  about  1300  km2,  and  the  area  of  WAD  is  about  2540  km2,  which  represent 
82%,  6%,  and  12%  of  the  computational  domain  within  Cl,  respectively. 

The  WAD  in  the  eastern  two  thirds  of  TA  depends  on  a  number  of  factors.  These  include  the 
bathymetry  in  the  western  third  of  TA  and,  in  the  eastern  two  thirds  of  TA,  the  bottom  depths 
and  slope,  the  channel  width  and  cross-sectional  shape,  and  the  bottom  drag.  There  is  a  significant 
amount  of  uncertainty  in  all  these  aspects  of  TA.  In  addition,  the  existence  of  small  drainage 
channels  (Ralston  and  Stacey  2007)  and  the  porosity  of  the  bottom  in  TA  are  not  accounted  for  in 
our  simulations  and  probably  affect  the  drainage  to  some  extent.  Because  of  these  uncertainties, 
trying  to  guess  the  proper  parameters  for  TA  is  a  process  that  could  go  on  indefinitely.  Hence,  we 
did  not  try  further  refinements  of  the  simulation  of  the  tidal  bore  in  TA  for  this  report. 

The  simulations  of  the  tidal  bore  in  TA  that  were  conducted  here  do  suggest  a  few  points.  The 
western  third  of  TA  must  be  sufficiently  deep  and  clear  of  obstructions  that  the  tidal  flow  entering 
the  eastern  two  thirds  of  TA  is  not  too  restricted  and  is  of  sufficient  volume.  The  bathymetry 
provided  to  us  by  Tal  Ezer  seems  to  be  too  shallow  and  restricted  in  the  western  third  of  TA.  Part 
of  the  problem  may  be  that  the  depths  are  referenced  to  something  lower  than  MSL,  and  need  to  be 
deepened.  The  ZP  bathymetry  as  received  was  referenced  to  MLLW.  In  adjusting  this  bathymetry 
to  MSL,  the  depths  in  the  entire  northeastern  part  of  Cl  were  deepened  significantly,  i.e.,  by  about 
5  nr.  The  use  of  this  bathymetry  resulted  in  a  strong,  rapid  response  of  the  tidal  flood  at  the 
entrance  to  the  narrow,  eastern  two  thirds  of  TA  with  respect  to  the  time  of  low  tide  at  Anchorage. 

Similarly,  the  bottom  depth  near  Beluga  Pt.  cannot  be  too  shallow.  This  would  seem  to  be 
indicated  by  the  fairly  rapid  appearance  of  the  tidal  bore  observed  at  Beluga  Pt.  after  the  time 
of  low  tide  at  Anchorage  (75  minutes).  If  the  bottom  elevation  near  Beluga  Pt.  is  too  high,  the 
tidal  flood  there  will  be  delayed  until  the  SSH  in  the  western  part  of  TA  rises  sufficiently  to  allow 
the  tide  to  enter  the  narrow,  eastern  two  thirds  of  TA.  The  effect  of  inertia  and  conservation  of 
momentum  can  push  the  tidal  height  near  Beluga  Pt  to  elevations  higher  than  those  in  the  western 
part  of  TA;  however,  the  simulations  we  conducted  suggest  that  this  effect  is  insufficient  to  match 
the  observed  time  of  the  tidal  bore  if  the  bottom  elevation  at  Beluga  Pt.  is  much  higher  than  about 
-3  m. 


If  the  tide  in  the  eastern  two  thirds  of  TA  is  to  drain,  then  there  must  be  a  fairly  constant 
upward  slope  from  west  to  east,  so  that  the  water  can  drain.  Observations  seem  to  indicate  that 
though  some  water  may  remain  behind,  most  of  the  water  in  the  eastern  two  thirds  of  TA  drains 
out  at  low  tide.  We  also  found  that,  for  most  of  the  water  to  drain  out  of  the  eastern  part  of  TA  at 
low  tide,  the  bottom  slope  must  be  sufficiently  steep.  To  get  some  time  period  of  complete  drainage 
at  Beluga  Pt.  on  the  lowest  tides  in  the  simulations  we  conducted  required  a  bottom  slope  of  about 
2  m  per  10  km. 
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3.8  Hurricane  Ike 

Hurricanes  are  associated  with  large  storm  surges  and  coastal  flooding.  These  processes  are 
important  to  capture  when  using  a  WAD  model.  In  this  section,  the  WAD  caused  by  Hurricane  Ike 
along  the  Texas  and  Louisiana  coasts  in  September  2008  is  simulated  and  the  results  are  compared 
with  observations. 

Ike  made  landfall  on  September  13,  2008  near  Galveston,  TX  as  a  category  2  hurricane  (Veer- 
arnony  et  al.  2014).  Ike  caused  significant  flooding  and  the  water  levels  during  the  storm  were 
observed  at  a  number  of  locations.  Hence,  Hurricane  Ike  provides  a  useful  test  case  for  the  verifi¬ 
cation  of  a  WAD  model. 

The  domain  and  bathymetry  used  for  the  simulations  of  Hurricane  Ike  are  shown  in  Fig.  18a. 
Figure  18b  shows  the  bathymetry  and  land  elevations  along  part  of  the  Texas  and  Louisiana  coasts 
up  to  a  height  of  10  m.  The  bathymetry  in  Fig.  18  is  the  same  as  that  used  by  Veeramony  et  al.  2014. 
The  bathymetry  was  derived  from  data  made  available  by  the  Southeastern  Universities  Research 
Association  (SURA)  Inundation  Testbed.  This  is  a  high-  (approx.  30-m)  resolution  dataset  that 
covers  much  of  the  northern  Gulf  of  Mexico.  In  the  deeper  water  not  covered  by  the  SURA  data, 
the  bathymetry  was  obtained  from  the  National  Geophysical  Data  Center  (NGDC)  Coastal  Relief 
Model,  Shuttle  Radar  Topography  Mission  (SRTM)  and  from  the  General  Bathymetric  Chart  of 
the  Ocean  (GEBCO). 

The  horizontal  grid  resolution  for  the  Ike  simulations  is  0.02°  (approximately  2  km).  The  grid 
size  is  504  x  348  points  in  the  x  (longitude)  and  y  (latitude)  directions,  respectively.  Grid  cells 
with  an  elevation  less  than  +10  m  were  taken  to  be  within  the  computational  domain  and,  hence, 
subject  to  WAD.  Grid  cells  with  an  elevation  >  +10  m  were  defined  to  be  land  points  outside  the 
computational  domain.  The  model  was  run  with  one  vertical  layer.  A  quadratic  form  of  the  bottom 
drag  was  used  (see  Section  2.5),  with  a  minimum  bottom  drag  coefficient  of  0.0025  and  bottom 
roughness  of  0.01  m.  A  third-order  upwind  scheme  was  used  for  momentum  advection.  The  time 
step  was  60  seconds.  Tidal  boundary  conditions  for  the  8  main  tidal  constituents  (Ki,  Or,  Pi,  Qi, 
K2,  M2,  N2,  S2)  from  the  OSU  1/45°  Gulf  of  Mexico  database  were  applied  at  the  open  boundaries. 

The  atmospheric  forcing  was  provided  by  Oceanweather  Inc.  (OWI)  on  a  longitude-latitude 
grid  covering  the  Gulf  of  Mexico  with  a  spatial  resolution  of  0.02°  and  a  temporal  resolution  of 
15  minutes,  and  included  atmospheric  pressure  and  wind  velocities.  These  fields  were  interpolated 
to  the  NCOM  grid,  and  the  winds  were  adjusted  for  land  effects  by  employing  a  directional  land- 
masking  scheme  (Veeramony  et  al.  2014;  Westerink  et  al.  2008).  Wind  stresses  were  computed  from 
the  wind  velocities  using  standard  bulk  formulas  with  the  wind-stress  drag  coefficient  Cd  computed 
as  Cd  =  0.00218  for  Ua  greater  than  1  m/s,  Cd  =  0.00062  +  0.00156 /Ua  for  Ua  between  1  and  3 
m/s,  Cd  =  0.00114  for  Ua  between  3  and  10  m/s,  Cy  =  0.00049  +  0.000065C/a  for  Ua  between  10 
and  26  m/s,  and  Cd  =  0.00216  for  Ua  greater  than  26  m/s,  where  Ua  is  the  wind  speed  in  m/s  (Xia 
et  al.  2008). 

The  Ike  simulations  were  started  from  rest  at  12Z  on  September  5,  2008.  This  is  well  before 
Ike  entered  the  Gulf  of  Mexico  on  September  9  and  allows  sufficient  time  for  the  model  to  spin-up 
before  Ike  enters  the  model  domain.  The  model  was  run  for  10  days  from  12Z  September  5  to  12Z 
September  15. 
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Fig.  18  -  -  Upper  plot  (a)  shows  model  domain  and  bathymetry  (m)  used  for  Hurricane  Ike  simulations.  Contours  are 
from  -3500  to  0  m.  Land  elevations  are  not  shown.  Lower  plot  (b)  shows  bathymetry  and  land  elevations  (m)  along 
a  portion  of  the  Texas  and  Louisiana  coasts.  Contours  are  from  -40  to  +10  m.  The  approximate  coastline  is  shown 
by  the  black  line.  The  locations  of  NOS  tide  stations  used  to  validate  the  SSH  predictions  are  indicated  by  red  dots. 
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Hurricane  Ike  Inundation 


Fig.  19  Hurricane  Ike  Inundation  Map.  (Obtained  from  Harris  County  Flood  Control  District  www.hcfcd.org) 


One  of  the  most  important  aspects  of  a  hurricane  that  people  in  the  path  of  the  storm  would  like 
to  know  is  where  and  when  flooding  will  occur  due  to  the  storm  surge.  The  estimated  inundation 
(depth  of  flooding)  caused  by  Ike,  as  determined  by  the  Harris  County  (Texas)  Flood  Control 
District,  is  shown  in  Fig.  19. 

To  qualitatively  compare  with  the  observed  inundation  map  in  Fig.  19,  the  maximum  inunda¬ 
tion  from  the  NCOM  simulation  is  shown  in  Fig.  20.  The  model  results  show  areas  of  maximum 
inundation  along  the  coast  from  Galveston  Bay,  TX  to  Mud  Lake,  LA.  The  inundation  from  the 
model  simulation  sometimes  reaches  further  inland  than  the  observed  inundation  in  Figure  19. 
However,  the  magnitude  of  the  observed  inundation  in  Fig.  19  appears  to  be  underestimated  by 
the  model  simulation  by  a  meter  or  more. 

During  Hurricane  Ike,  several  National  Ocean  Service  (NOS)  tide  stations  along  the  Texas- 
Louisiana  coast  recorded  the  water  levels  as  the  storm  made  landfall.  Figure  18b  shows  the  locations 
of  six  NOS  stations  used  to  compare  the  observed  water  levels  with  the  model  results. 

Figure  21  contains  plots  comparing  the  NOS  observed  water  levels  (black  lines)  with  the  NCOM 
simulated  water  levels  (blue  lines).  In  this  figure,  it  can  be  seen  that,  up  to  September  9,  the  ob¬ 
served  and  simulated  water  levels  are  relatively  close,  with  the  model  under-estimating  the  observed 
water  levels  by  about  0.2  m.  After  September  9,  as  Ike  begins  to  move  closer  to  shore  and  affect 
the  water  levels  near  shore,  the  model  water  levels  rise,  but  under-estimate  the  observed  values  by 
about  a  meter.  One  possible  explanation  for  the  under-estimation  of  the  water  levels  is  that  the 
effect  of  waves  on  the  water  levels  is  not  accounted  for.  Veeramony  et  al.  (2014)  found  that  the 
inclusion  of  wave  processes  increased  the  modeled  water  levels  by  an  average  of  about  0.8  m  (see 
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Fig.  20  -  Maximum  water  inundations  (m)  from  NCOM  simulation. 

their  Fig.  2.1-15).  If  the  inclusion  of  wave  processes  in  NCOM  resulted  in  a  similar  increase  in  the 
mean  water  levels,  the  comparison  with  the  observed  water  levels  would  be  significantly  improved. 
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A  baroclinic,  multi-layer  simulation  of  Hurricane  Ike  was  run  using  a  vertical  grid  with  a  total 
of  40  layers,  with  19  sigma  layers  between  the  surface  and  the  bottom  down  to  a  depth  of  -115  m 
and  21  fixed-depth  levels  between  -115  m  and  -3500  m.  Hence,  the  vertical  grid  in  water  shallower 
than  -115  m,  including  the  WAD  areas,  was  19  sigma  layers.  The  maximum  static  surface  layer 
thickness  in  deep  water  was  1  m  and  the  grid  was  smoothly  stretched  in  the  vertical,  with  each 
layer  being  17%  thicker  than  the  layer  above.  The  horizontal  grid  and  bathymetry  were  the  same 
as  used  for  the  barotropic  simulation  (Fig.  18).  The  SSH,  velocity,  temperature,  and  salinity  were 
initialized  from  operational  Global  NCOM  fields  obtained  from  the  Naval  Oceanographic  Office 
(NAVO).  The  lateral  boundary  conditions  consisted  of  tidal  forcing  as  used  for  the  barotropic 
simulation  plus  daily  values  of  SSH,  velocity,  temperature,  and  salinity  from  Global  NCOM.  The 
surface  atmospheric  forcing  was  the  surface  pressure  and  wind  stresses  from  Oceanweather  Inc.  as 
used  for  the  barotropic  simulation.  The  timestep  used  was  60  s,  the  same  as  for  the  barotropic 
simulation.  This  baroclinic,  multi-layer,  WAD  simulation  of  Hurricane  Ike  ran  smoothly.  Figure  22 
shows  the  model-predicted  SST  at  00Z  September  14  after  Ike  had  gone  ashore  near  Galveston,  TX. 
The  SST  shows  surface  cooling  of  up  to  5°C  along  the  track  of  the  hurricane  caused  by  upwelling 
and  vertical  mixing. 

4.  LIMITATIONS  OF  WAD  IN  NCOM 
4.1  Time-Step  Limitations 

The  WAD  results  in  NCOM  were  not  found  to  be  very  sensitive  to  the  time  step,  despite 
that  fact  that  a  relatively  large  time  step  is  used  for  the  implicit  solution  of  the  barotropic  mode 
in  NCOM.  This  is  probably  because  the  WAD  scheme  in  NCOM  always  maintains  a  minimum 
specified  water  depth  Dmin  at  all  the  grid  cells  within  the  computational  domain.  This  avoids 
having  to  reduce  the  time  step  to  prevent  the  water  depth  D  from  dropping  from  above  Drnin  to 
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Fig.  21  —  Observed  NOS  (black)  and  model-simulated  (blue)  water  levels  for  Hurricane  Ike.  Note  that  the  NOS  gage 
at  Galveston  North  Jetty  stopped  reporting  values  on  September  13.  The  flat  blue  line  at  Eagle  Point  indicates  that 
this  area  was  dry  in  the  model  until  September  13. 
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Fig.  22  —  Model-simulated  SST  for  Hurricane  Ike  at  00Z  September  14,  2008  in  °C. 
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below  the  bottom  during  a  given  time  step.  However,  the  use  of  a  larger  time  step  may  increase  the 
number  of  times  the  barotropic  solver  must  be  rerun  during  a  single  time  step  if  there  are  extensive 
WAD  areas.  This  is  because,  for  a  larger  time  step,  there  tend  to  be  larger  changes  in  the  water 
depths  for  the  WAD  scheme  to  deal  with. 

However,  the  time  step  must  be  small  enough  to  avoid  exceeding  the  CFL  constraints  for 
advection.  The  occurrence  of  small  water  depths  due  to  specifying  a  small  value  of  Dmin  and  the 
presence  of  relatively  steep  bottom  slopes  in  WAD  areas  can  both  result  in  increased  velocities  that 
may  require  reduction  of  the  time  step  to  avoid  exceeding  the  CFL  limits  for  advection. 


4.2  Bathymetry  Limitations 

There  are  no  explicit  limitations  on  the  bathymetry  for  WAD  in  NCOM.  However,  as  noted 
previously,  the  presence  of  steep  bottom  slopes  can  generate  larger  vertical  velocities  that  may 
require  a  reduction  of  the  time  step  to  avoid  violating  the  CFL  constraint  for  vertical  advection. 
Hence,  having  smooth  changes  in  the  bottom  depth  in  WAD  areas  may  allow  the  use  of  a  larger 
time  step  and,  consequently,  reduce  the  run  time. 

4.3  WAD  near  Open  Boundaries 

At  the  present  time,  WAD  cannot  be  handled  at  the  open  boundaries.  Hence,  a  maximum 
allowable  static  bottom  depth  H  must  be  specified  at  grid  cells  at  open  boundaries  within  the 
computational  domain  so  WAD  will  not  occur  there.  WAD  at  the  open  boundaries  is  currently 
being  investigated  in  order  to  try  to  remove  this  restriction. 

4.4  Effect  of  WAD  on  Run  Time 

The  use  of  WAD  in  NCOM  can  increase  the  run  time,  both  by  increasing  the  time  required 
during  a  single  time  step  (due  to  the  need  to  recompute  the  barotropic  solver  in  NCOM  a  number 
of  times  during  a  single  time  step  when  there  is  WAD),  and  by  increasing  the  number  of  time  steps 
that  are  required  (due  to  a  need  to  decrease  the  time  step  to  avoid  violating  the  CFL  constraints 
for  advection). 

The  increase  in  the  run  time  for  “occasional”  WAD,  such  as  was  used  in  the  ChB  simulations 
described  in  this  report,  tends  to  be  small.  As  noted  previously,  “occasional”  WAD  refers  to  WAD 
that  occurs  only  occasionally  in  a  domain  in  which  some  maximum  allowable  static  bottom  depth 
Hmaxi  with  Hmax  <  0  (e.g.,  H  <  Hmax  =  —2  nr),  has  been  prescribed  for  all  the  grid  cells  within 
the  computational  domain. 

For  “extensive”  WAD,  where  there  are  extensive  WAD  areas  and  WAD  is  almost  always  oc¬ 
curring  somewhere,  such  as  in  the  SFB  and  Cl  simulations  described  in  this  report,  the  increase  in 
the  run  time  will  tend  to  be  more  significant  (due  to  the  need  to  recompute  the  barotropic  solver 
in  NCOM  a  number  of  times  during  every,  or  almost  every,  time  step).  And  as  the  grid  resolution 
is  increased  in  areas  where  there  is  “extensive”  WAD,  the  increase  in  the  run  time  needed  for  a 
single  time  step  will  tend  to  increase  further. 


45 


Martin  et  al. 


5.  SUMMARY 

This  report  discusses  the  implementation  and  testing  of  wetting  and  drying  (WAD)  in  the  Navy 
Coastal  Ocean  Model  (NCOM).  NCOM  is  run  as  a  stand-alone  ocean  model  and  also  as  part  of  the 
Coupled  Ocean/ Atmosphere  Mesoscale  Prediction  System  (COAMPS),  which  provides  for  one-  or 
two-way  coupling  among  atmosphere,  ocean,  and  wave  models. 

The  implementation  of  WAD  in  NCOM  was  complicated  by  the  fact  that  NCOM  uses  an 
implicit  numerical  scheme  to  update  the  barotropic  mode,  i.e. ,  to  update  the  SSH  and  the  depth- 
integrated  transports.  With  the  use  of  an  implicit  scheme  to  update  the  barotropic  mode,  the 
same,  relatively  large,  timestep  is  used  for  the  update  of  the  barotropic  mode  that  is  used  for  the 
rest  of  the  ocean  model.  This  is  possible  because  an  implicit  treatment  of  the  barotropic  mode  is 
not  limited  by  the  speed  of  surface  gravity  waves,  and  so  a  much  larger  timestep  can  be  used  than 
if  the  barotropic  mode  were  updated  explicitly. 

The  WAD  scheme  is  primarily  implemented  within  NCOM’s  solution  of  the  barotropic  mode. 
After  the  new  SSH  is  computed  (using  an  iterative  solver),  the  new  water  depth  in  each  grid  cell 
is  inspected  and,  for  grid  cells  whose  depths  have  fallen  below  a  prescribed  minimum  Dmin ,  the 
volume  fluxes  at  cell  faces  that  are  directed  out  of  the  drying  grid  cells  are  set  to  zero  by  setting 
the  solver  coefficient  for  those  grid-cell  faces  to  zero,  and  then  the  free-surface  solver  is  rerun.  This 
procedure  is  repeated  until  convergence  occurs,  i.e.,  until  no  grid-cell  depths  drop  below  Drnin  or 
the  maximum  difference  in  the  SSH  from  the  calculation  of  the  SSH  on  the  previous  run  of  the 
solver  falls  below  a  small  prescribed  value.  Hence,  with  this  scheme,  the  water  depth  at  grid  cells 
within  the  computational  domain  is  not  allowed  to  fall  much  below  the  minimum  specified  depth 

Dmin- 


The  WAD  in  NCOM  was  tested  by  running  simulations  of  (a)  idealized  experiments  that  have 
expected  or  analytical  solutions  that  can  be  compared  against,  (b)  laboratory  experiments  that 
have  observed  results  that  can  be  compared  against,  (c)  several  coastal  regions  that  have  notable 
WAD  areas,  i.e.,  San  Francisco  Bay  (SFB),  Chesapeake  Bay  (ChB),  and  Cook  Inlet  in  Alaska,  and 
(d)  Hurricane  Ike,  which  caused  extensive  flooding  along  the  Texas  and  Louisiana  coasts  in  2008. 

The  first  of  the  idealized  experiments  consisted  of  a  test  of  NCOM’s  horizontal  symmetry  when 
there  are  WAD  areas  within  the  domain  and  WAD  is  occurring.  This  is  a  test  of  NCOM’s  ability  to 
exactly  maintain  a  horizontally  symmetric  solution  for  a  horizontally  symmetric  problem.  This  is 
one  of  the  first  tests  conducted  with  NCOM  after  extensive  changes  have  been  made,  since  this  type 
of  test  can  detect  many  types  of  programming  errors.  NCOM  maintained  a  perfectly  symmetric 
solution  for  this  test. 

The  second  of  the  idealized  tests  consisted  of  the  propagation  of  a  planar  surface  wave  rotat¬ 
ing  around  within  a  parabolic-shaped  basin.  Wetting  occurs  at  the  leading  edge  of  the  wave  as 
it  propagates  around  the  basin,  and  drying  occurs  at  the  trailing  edge.  The  NCOM  simulation 
was  compared  with  the  analytical  solution.  The  NCOM  solution  agreed  with  the  analytical  solu¬ 
tion,  maintained  the  correct  period  of  rotation  of  the  wave  around  the  basin,  and  maintained  the 
amplitude  of  the  wave  for  several  revolutions  with  very  little  damping. 

The  first  of  the  laboratory  experiments  consisted  of  a  dam-break  flow  over  a  triangular  bump. 
Water  in  a  tank  is  released  at  the  beginning  of  the  experiment  and  flows  down  a  channel  in  which 
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there  is  a  triangular-shaped  bump.  Some  of  the  water  flows  over  the  bump,  and  a  surface  wave  is 
generated  at  the  bump  that  propagates  back  towards  the  tank,  reflects  from  the  back  wall  of  the 
tank,  and  propagates  back  towards  the  bump.  This  sequence  is  repeated  several  times  until  the 
surface  wave  propagating  between  the  bump  and  the  back  of  the  tank  is  dissipated.  The  water 
height  from  the  simulation  was  compared  with  the  observed  water  height  at  seven  different  locations 
and  the  agreement  was  fairly  good  at  all  seven  locations. 

The  second  of  the  laboratory  experiments  consisted  of  a  dam-break  flow  down  a  channel  and 
around  a  90°  bend  in  the  channel.  At  the  bend,  some  of  the  water  makes  its  way  around  the  bend 
and  some  is  reflected  back  towards  the  tank.  The  observed  water  height  was  measured  during  the 
experiment  at  six  different  locations,  and  the  agreement  between  the  simulated  water  height  and 
the  measured  water  height  was  fairly  good  at  all  of  the  locations. 

The  SFB  simulations  demonstrate  that  the  WAD  in  NCOM  can  flood  and  drain  extensive  areas 
when  used  in  different  model  configurations,  including  different  grid  resolutions  (500,  200,  and  100 
m),  in  baroclinic  or  barotropic  mode,  with  different  numbers  of  layers  (1-40),  and  with  different 
numerical  options. 

The  ChB  case  was  run  to  illustrate  the  case  of  occasional  WAD,  where  WAD  occurs  only 
occasionally  due  to  drying  of  grid  cells  caused  by  very  low  tides  and/or  strong,  offshore  winds. 
Without  a  WAD  capability,  the  drying  out  of  grid  cells  caused  NCOM  to  crash.  With  the  WAD 
capability,  NCOM  continues  to  run  smoothly  when  grid  cells  dry  out.  With  just  occasional  WAD, 
there  is  no  significant  increase  in  the  NCOM  run  time. 

The  Cook  Inlet  case  demonstrates  that  the  WAD  scheme  performs  well  in  a  region  of  very  large 
(greater  than  10  m)  tidal  range  and  extensive  WAD  areas.  Cook  Inlet  has  large  areas  of  mud  flats 
that  are  exposed  at  low  tide.  We  focused  on  the  tidal  bore  and  flood  that  occurs  in  Turnagain 
Arm  near  the  head  of  Cook  Inlet.  The  eastern  two-thirds  of  Turnagain  Arm,  which  has  a  length  of 
about  50  km,  almost  fully  drains  at  low  tide.  The  tidal  bore  that  occurs  when  the  tidal  flood  enters 
Turnagain  Arm  is  a  popular  attraction;  hence,  the  approximate  arrival  times  of  the  bore  at  various 
locations  along  Turnagain  Arm  are  well  known.  Initially,  our  simulated  flood  of  Turnagain  Arm 
occurred  much  later  than  observed.  However,  with  improved  bathymetry,  we  were  able  to  achieve 
times  that  were  much  closer  to  the  observed  times. 

The  Hurricane  Ike  case  demonstrates  that  the  WAD  performs  well  in  high  wind  conditions  in 
which  the  WAD  is  caused  by  large  storm  surge.  However,  the  NCOM  simulation  underestimated 
the  observed  storm  surge  by  about  a  meter.  This  underestimate  may  be  due,  at  least  in  part,  to 
the  neglect  of  wave  effects  on  the  storm  surge. 

The  pros  and  cons  of  the  WAD  scheme  implemented  in  NCOM  are  discussed.  The  WAD  scheme 
has  the  main  advantages  that  it  does  not  require  any  special  modification  of  the  bathymetry  and  is 
fairly  robust.  The  main  disadvantages  of  the  WAD  scheme  are  that  (a)  areas  subject  to  WAD  must 
be  defined  to  be  within  the  computational  domain  at  the  start  of  the  simulation,  (b)  a  minimum 
thickness  of  fluid  must  exist  at  all  times  within  the  WAD  areas,  i.e.,  the  WAD  areas  can  never  be 
totally  dry,  and  (c)  additional  calculations  are  required  for  the  WAD  and  the  timestep  may  need 
to  be  decreased,  both  of  which  increase  the  NCOM  run  time. 
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